Coral Performance Following Reciprocal Transplant

Setup

Load Packages

##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("cowplot")) install.packages("cowplot")
if (!require("lubridate")) install.packages("lubridate")

Loading required package: lubridate

Attaching package: ‘lubridate’

The following object is masked from ‘package:cowplot’:

    stamp

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

if (!require("corrplot")) install.packages("corrplot")
if (!require("lme4")) install.packages("lme4")
if (!require("vegan")) install.packages("vegan")
if (!require("DHARMa")) install.packages("DHARMa")
if (!require("effectsize")) install.packages("effectsize")
if (!require("emmeans")) install.packages("emmeans")
if (!require("lmerTest")) install.packages("lmerTest")
if (!require("tidyr")) install.packages("tidyr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("plotrix")) install.packages("plotrix")

Loading required package: plotrix
Warning: package ‘plotrix’ was built under R version 4.3.2

if (!require("Rmisc")) install.packages("Rmisc")
if (!require("ggpubr")) install.packages("ggpubr")

Loading required package: ggpubr

Attaching package: ‘ggpubr’

The following object is masked from ‘package:plyr’:

    mutate

The following object is masked from ‘package:cowplot’:

    get_legend

##Load Packages
library(ggplot2) #Required for ggplots
library(cowplot) #Required for plotting panel figures
library(lubridate) #Required for date and time formatting
library(corrplot) #Required for correlation plot
library(lme4) #Required for mixed effects modeling
library(vegan) #Required for multivariate analysis PERM
library(DHARMa) #Required to check residuals of mixed effects models
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons 
library(lmerTest) #Required for p values with lme4 model summaries
library(tidyr) #Required for pivot_wider function
library(dplyr) #Required for group_by function
library(plotrix) #Required for Standard error function
library(Rmisc) #Required for SummarySE function
library(ggpubr) #Required for stat brackets

Graphing Parameters

#Note: Run "Graphing Parameters" section from 01_ExperimentalSetup.Rmd file

Sample Data and Metadata

Load and Organize Data

##Growth
Length<-read.csv("Data/LengthFull.csv", header=TRUE)
Dates<-read.csv("Data/BonaireDates.csv", header=TRUE)

##Set date variables
Dates$Date<-as.Date(Dates$Date, format='%m/%d/%Y')

##Set factor variables
Length$Site<-factor(Length$Site, levels=c("KL", "SS"))
Length$Genotype<-factor(Length$Genotype, levels=c("AC8", "AC10", "AC12"))
Length$Orig<-factor(Length$Orig, levels=c("N", "T"))
Length$Origin<-factor(Length$Origin, levels=c("Native", "Transplant"))

##Add a Sample Set Variable
Length$Set<-paste(Length$Site, Length$Genotype, Length$Orig, sep=".")
Length$Set<-factor(Length$Set)

##Add Site.Orig variable
Length$Site.Orig<-paste(Length$Site, Length$Orig, sep=".")
Length$Site.Orig<-factor(Length$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))


##Bleaching Metrics
#Note: Physiological metrics calculated in 02_PhysiologyMetrics.R file
Thermal<-read.csv("Outputs/ThermData.csv", header=TRUE)

##Set factor variables
Thermal$TimeP<-factor(Thermal$TimeP, levels=c("IN", "TP1", "TP2", "TP3", "TP4"))
Thermal$Site<-factor(Thermal$Site, levels=c("KL", "SS"))
Thermal$Genotype<-factor(Thermal$Genotype, levels=c("AC8", "AC10", "AC12"))
Thermal$Orig<-factor(Thermal$Orig, levels=c("N", "T"))
Thermal$Origin<-factor(Thermal$Origin, levels=c("Native", "Transplant"))

Growth

Growth Rate TP1 to TP2

Calculate Days

##Subset Dates for TP 1 and TP 2 Growth
Growth.Dates_T1.2<-subset(Dates, Task=="Growth" & TimeP=="TP1" |
                       Task=="Growth" & TimeP=="TP2" )

##Convert to Wide format
Growth.Dates_T1.2<-Growth.Dates_T1.2 %>%  pivot_wider(names_from = TimeP, values_from = Date)

##Calculate Number of Days
#TP 1 to TP2
Growth.Dates_T1.2$TP1.2_days<-time_length(Growth.Dates_T1.2$TP2-Growth.Dates_T1.2$TP1, unit="days")

Calculate Total Linear Extension

##Merge Length Data with Growth Dates
#Merges by Site column
#Adds TP1.2_days column
Growth_T1.2<-merge(Length, Growth.Dates_T1.2[,c(-2)], all.x=TRUE)

##Calculate Linear Extension (cm)
Growth_T1.2$Ext_cm<-Growth_T1.2$TL.TP2_cm-Growth_T1.2$TL.TP1_cm

##Calculate Growth_T1.2 Rate (cm/day)
Growth_T1.2$TLE_cm.day<-Growth_T1.2$Ext_cm/Growth_T1.2$TP1.2_days

Data QC

##Remove Corals not Measured in T1 to T2
Growth_T1.2<-Growth_T1.2 %>% drop_na(TLE_cm.day)

##Sample Size per Set
Growth_T1.2 %>% 
  dplyr::group_by(Set) %>%
  dplyr::summarize(SampleSize = length(Set))

##Check for Outliers
ggplot(Growth_T1.2, aes(x=Set, y=TLE_cm.day)) + 
  geom_boxplot(alpha=0.5, shape=2, outlier.shape = NA)+
  geom_jitter(shape=16, position=position_jitter(0.1))+
  theme(axis.text.x = element_text(angle = 90))

Potential outliers have been re-measured and therefore will be retained. n=10 per Site, Genotype, Origin group

Compare Growth Rate

#Check normality
hist(Growth_T1.2$TLE_cm.day)

shapiro.test(Growth_T1.2$TLE_cm.day)

    Shapiro-Wilk normality test

data:  Growth_T1.2$TLE_cm.day
W = 0.92803, p-value = 7.262e-06

#Not normal

hist(log(Growth_T1.2$TLE_cm.day+1))

shapiro.test(log(Growth_T1.2$TLE_cm.day+1))

    Shapiro-Wilk normality test

data:  log(Growth_T1.2$TLE_cm.day + 1)
W = 0.94335, p-value = 7.248e-05

##Still not normal

##Compare generalized linear mixed effects models
TLE_T1.2.glmr.gaus<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=gaussian(link="identity"))

Warning: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly

TLE_T1.2.glmr.gaus.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=gaussian(link="log"))
TLE_T1.2.glmr.gaus.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=gaussian(link="inverse"))

TLE_T1.2.glmr.gam.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=Gamma(link="log"))
TLE_T1.2.glmr.gam.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=Gamma(link="inverse"))

TLE_T1.2.glmr.inv<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=inverse.gaussian(link="1/mu^2"))

boundary (singular) fit: see help('isSingular')

TLE_T1.2.glmr.inv.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=inverse.gaussian(link="log"))
TLE_T1.2.glmr.inv.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T1.2, family=inverse.gaussian(link="inverse"))

AIC(TLE_T1.2.glmr.gaus, TLE_T1.2.glmr.gaus.l, TLE_T1.2.glmr.gaus.i, TLE_T1.2.glmr.gam.l, TLE_T1.2.glmr.gam.i, TLE_T1.2.glmr.inv, TLE_T1.2.glmr.inv.l, TLE_T1.2.glmr.inv.i)

Gamma distribution with log-link has the lowest AIC.

##Check residuals
TLE_T1.2.glmr.gam.l_res <- simulateResiduals(fittedModel = TLE_T1.2.glmr.gam.l, plot = F)
plot(TLE_T1.2.glmr.gam.l_res)

Compare residuals across models

##Check residuals
TLE_T1.2.glmr.gaus.l_res <- simulateResiduals(fittedModel = TLE_T1.2.glmr.gaus.l, plot = F)
plot(TLE_T1.2.glmr.gaus.l_res)

Compare with log+1 transformed lmer model

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
TLE_T1.2.lme<-lmer(log(TLE_cm.day+1)~Origin*Site+(1|Genotype), data=Growth_T1.2)

##Check residuals
TLE_T1.2.lme_res <- simulateResiduals(fittedModel = TLE_T1.2.lme, plot = F)
plot(TLE_T1.2.lme_res)

Better residuals with the log+1 transformed lmer model.

##Model results
summary(TLE_T1.2.lme)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin * Site + (1 | Genotype)
   Data: Growth_T1.2

REML criterion at convergence: -233.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0026 -0.5991 -0.1118  0.5959  2.4282 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002965 0.05445 
 Residual             0.006632 0.08144 
Number of obs: 120, groups:  Genotype, 3

Fixed effects:
                          Estimate Std. Error         df t value Pr(>|t|)   
(Intercept)               0.125179   0.034778   2.684719   3.599  0.04400 * 
OriginTransplant         -0.007538   0.021027 114.000000  -0.359  0.72063   
SiteSS                    0.062096   0.021027 114.000000   2.953  0.00382 **
OriginTransplant:SiteSS  -0.035167   0.029737 114.000000  -1.183  0.23943   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.302              
SiteSS      -0.302  0.500       
OrgnTrn:SSS  0.214 -0.707 -0.707

eta_squared(TLE_T1.2.lme)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.02 | [0.00, 1.00]
Site        |           0.07 | [0.01, 1.00]
Origin:Site |           0.01 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
TLE_T1.2.res<-data.frame(summary(TLE_T1.2.lme)$coefficients[-1,])
TLE_T1.2.res$Predictor<-c("Origin", "Site", "Origin x Site")
TLE_T1.2.res$EtaSq<-c(eta_squared(TLE_T1.2.lme)$Eta2)
TLE_T1.2.res$Response<-rep("Growth", nrow(TLE_T1.2.res))
TLE_T1.2.res$Timepoint<-rep("TP1_2", nrow(TLE_T1.2.res))

Effect size of Origin for each Site

##KL
TLE_T1.2.lme_KL<-lmer(log(TLE_cm.day+1)~Origin+(1|Genotype), data=Growth_T1.2[which(Growth_T1.2$Site=="KL"),])
summary(TLE_T1.2.lme_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin + (1 | Genotype)
   Data: Growth_T1.2[which(Growth_T1.2$Site == "KL"), ]

REML criterion at convergence: -115.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7586 -0.6776 -0.1453  0.5574  2.2922 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.001699 0.04122 
 Residual             0.006724 0.08200 
Number of obs: 60, groups:  Genotype, 3

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)       0.125179   0.028116  2.712715   4.452    0.026 *
OriginTransplant -0.007538   0.021173 56.000006  -0.356    0.723  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.377

eta_squared(TLE_T1.2.lme_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |       2.26e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
TLE_T1.2.lme_SS<-lmer(log(TLE_cm.day+1)~Origin+(1|Genotype), data=Growth_T1.2[which(Growth_T1.2$Site=="SS"),])
summary(TLE_T1.2.lme_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin + (1 | Genotype)
   Data: Growth_T1.2[which(Growth_T1.2$Site == "SS"), ]

REML criterion at convergence: -114.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.05792 -0.57732 -0.09493  0.72022  2.38351 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.004144 0.06437 
 Residual             0.006600 0.08124 
Number of obs: 60, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.18727    0.04002  2.30548   4.680   0.0323 *
OriginTransplant -0.04271    0.02098 56.00000  -2.036   0.0465 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.262

eta_squared(TLE_T1.2.lme_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
TLE_T1.2.site.res<-data.frame(rbind(summary(TLE_T1.2.lme_KL)$coefficients[-1,],
                                     summary(TLE_T1.2.lme_SS)$coefficients[-1,]))
TLE_T1.2.site.res$Predictor<-c("KL Origin", "SS Origin")
TLE_T1.2.site.res$EtaSq<-c(eta_squared(TLE_T1.2.lme_KL)$Eta2, eta_squared(TLE_T1.2.lme_SS)$Eta2)
TLE_T1.2.site.res$Response<-rep("Growth", nrow(TLE_T1.2.site.res))
TLE_T1.2.site.res$Timepoint<-rep("TP1_2", nrow(TLE_T1.2.site.res))

##Combine results 
TLE_T1.2.res<-rbind(TLE_T1.2.res, TLE_T1.2.site.res)

##Summary statistics by Site and Origin
TP1.2_Growth.sum<-summarySE(Growth_T1.2, measurevar="TLE_cm.day", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Growth Rate across Treatments
TP1.2_Growth.plot<-ggplot(TP1.2_Growth.sum, aes(x=Site, y=TLE_cm.day, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=TLE_cm.day-se, ymax=TLE_cm.day+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y=expression(paste('Growth Rate (cm day'^-1*")")), colour="Origin")+
  ylim(0, 0.25)+
  annotate("text", x=2, y=0.24, label="*", size=sig.sz, fontface="bold")+
  geom_bracket(xmin=1, xmax=2, y.position=0.01, label.size=levels.sz, label="**", inherit.aes = FALSE); TP1.2_Growth.plot

Average Growth Rates by Genotype

##Summary statistics by Genotype and Origin
TP1.2_Growth.sum.SS<-summarySE(subset(Growth_T1.2, Site=="SS"), measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin"), na.rm=TRUE)

TP1.2_Growth.sum.SS

##Calculate Percent Difference between Native and Transplant by Genotype
TP1.2_Growth.SS<-TP1.2_Growth.sum.SS[which(TP1.2_Growth.sum.SS$Origin=="Native"), c(1,4)]
names(TP1.2_Growth.SS)[2]<-paste(names(TP1.2_Growth.SS)[2], "N", sep="_")
TP1.2_Growth.SS$TLE_cm.day_T<-TP1.2_Growth.sum.SS$TLE_cm.day[which(TP1.2_Growth.sum.SS$Origin=="Transplant")]
TP1.2_Growth.SS$Perc.Inc<-((TP1.2_Growth.SS$TLE_cm.day_N-TP1.2_Growth.SS$TLE_cm.day_T)/TP1.2_Growth.SS$TLE_cm.day_T)*100

TP1.2_Growth.SS

##Mean and Standard Error
mean(TP1.2_Growth.SS$Perc.Inc)

[1] 47.19841

std.error(TP1.2_Growth.SS$Perc.Inc)

[1] 25.89777

Growth Rate TP3 to TP4

Calculate Days

##Subset Dates for TP 3 and TP 4 Growth
Growth.Dates_T3.4<-subset(Dates, Task=="Growth" & TimeP=="TP3" |
                       Task=="Growth" & TimeP=="TP4" )

##Convert to Wide format
Growth.Dates_T3.4<-Growth.Dates_T3.4 %>%  pivot_wider(names_from = TimeP, values_from = Date)

##Calculate Number of Days
#TP 3 to TP4
Growth.Dates_T3.4$TP3.4_days<-time_length(Growth.Dates_T3.4$TP4-Growth.Dates_T3.4$TP3, unit="days")

Calculate Total Linear Extension

##Merge Length Data with Growth Dates
#Merges by Site column
#Adds TP3.4_days column
Growth_T3.4<-merge(Length, Growth.Dates_T3.4[,c(-2)], all.x=TRUE)

##Calculate Linear Extension (cm)
Growth_T3.4$Ext_cm<-Growth_T3.4$TL.TP4_cm-Growth_T3.4$TL.TP3_cm

##Calculate Growth_T3.4 Rate (cm/day)
Growth_T3.4$TLE_cm.day<-Growth_T3.4$Ext_cm/Growth_T3.4$TP3.4_days

Data QC

##Remove Corals not Measured in T1 to T2
Growth_T3.4<-Growth_T3.4 %>% drop_na(TLE_cm.day)

##Sample Size per Set
Growth_T3.4 %>% 
  dplyr::group_by(Set) %>%
  dplyr::summarize(SampleSize = length(Set))

##Check for Outliers
ggplot(Growth_T3.4, aes(x=Set, y=TLE_cm.day)) + 
  geom_boxplot(alpha=0.5, shape=2, outlier.shape = NA)+
  geom_jitter(shape=16, position=position_jitter(0.1))+
  theme(axis.text.x = element_text(angle = 90))

Compare Growth Rate

#Check normality
hist(Growth_T3.4$TLE_cm.day)

shapiro.test(Growth_T3.4$TLE_cm.day)

    Shapiro-Wilk normality test

data:  Growth_T3.4$TLE_cm.day
W = 0.92767, p-value = 6.904e-06

#Not normal

hist(log(Growth_T3.4$TLE_cm.day+1))

shapiro.test(log(Growth_T3.4$TLE_cm.day+1))

    Shapiro-Wilk normality test

data:  log(Growth_T3.4$TLE_cm.day + 1)
W = 0.94719, p-value = 0.0001351

##Still not normal

##Compare generalized linear mixed effects models
TLE_T3.4.glmr.gaus<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=gaussian(link="identity"))

Warning: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly

TLE_T3.4.glmr.gaus.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=gaussian(link="log"))
TLE_T3.4.glmr.gaus.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=gaussian(link="inverse"))

TLE_T3.4.glmr.gam.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=Gamma(link="log"))
TLE_T3.4.glmr.gam.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=Gamma(link="inverse"))

TLE_T3.4.glmr.inv<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=inverse.gaussian(link="1/mu^2"))
TLE_T3.4.glmr.inv.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=inverse.gaussian(link="log"))
TLE_T3.4.glmr.inv.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_T3.4, family=inverse.gaussian(link="inverse"))

AIC(TLE_T3.4.glmr.gaus, TLE_T3.4.glmr.gaus.l, TLE_T3.4.glmr.gaus.i, TLE_T3.4.glmr.gam.l, TLE_T3.4.glmr.gam.i, TLE_T3.4.glmr.inv, TLE_T3.4.glmr.inv.l, TLE_T3.4.glmr.inv.i)

Gamma distribution with log-link has the lowest AIC.

##Check residuals
TLE_T3.4.glmr.gam.l_res <- simulateResiduals(fittedModel = TLE_T3.4.glmr.gam.l, plot = F)
plot(TLE_T3.4.glmr.gam.l_res)

Compare residuals across models

##Check residuals
TLE_T3.4.glmr.gaus.l_res <- simulateResiduals(fittedModel = TLE_T3.4.glmr.gaus.l, plot = F)
plot(TLE_T3.4.glmr.gaus.l_res)

Compare with log+1 transformed lmer model

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
TLE_T3.4.lme<-lmer(log(TLE_cm.day+1)~Origin*Site+(1|Genotype), data=Growth_T3.4)

##Check residuals
TLE_T3.4.lme_res <- simulateResiduals(fittedModel = TLE_T3.4.lme, plot = F)
plot(TLE_T3.4.lme_res)

Better residuals with the log+1 transformed lmer model.

##Model results
summary(TLE_T3.4.lme)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin * Site + (1 | Genotype)
   Data: Growth_T3.4

REML criterion at convergence: -341.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.27432 -0.60263 -0.01933  0.56722  2.45603 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002008 0.04481 
 Residual             0.002574 0.05074 
Number of obs: 120, groups:  Genotype, 3

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)               0.16891    0.02748   2.38972   6.147 0.016262 *  
OriginTransplant         -0.02845    0.01310 114.00000  -2.172 0.031936 *  
SiteSS                   -0.05868    0.01310 114.00000  -4.479 1.79e-05 ***
OriginTransplant:SiteSS   0.06284    0.01853 114.00000   3.392 0.000954 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.238              
SiteSS      -0.238  0.500       
OrgnTrn:SSS  0.169 -0.707 -0.707

eta_squared(TLE_T3.4.lme)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       9.01e-04 | [0.00, 1.00]
Site        |           0.07 | [0.01, 1.00]
Origin:Site |           0.09 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
TLE_T3.4.res<-data.frame(summary(TLE_T3.4.lme)$coefficients[-1,])
TLE_T3.4.res$Predictor<-c("Origin", "Site", "Origin x Site")
TLE_T3.4.res$EtaSq<-c(eta_squared(TLE_T3.4.lme)$Eta2)
TLE_T3.4.res$Response<-rep("Growth", nrow(TLE_T3.4.res))
TLE_T3.4.res$Timepoint<-rep("TP3_4", nrow(TLE_T3.4.res))

Effect size of Origin for each Site

##KL
TLE_T3.4.lme_KL<-lmer(log(TLE_cm.day+1)~Origin+(1|Genotype), data=Growth_T3.4[which(Growth_T3.4$Site=="KL"),])
summary(TLE_T3.4.lme_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin + (1 | Genotype)
   Data: Growth_T3.4[which(Growth_T3.4$Site == "KL"), ]

REML criterion at convergence: -155.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.13557 -0.60999 -0.09472  0.50939  2.13461 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002615 0.05113 
 Residual             0.003233 0.05686 
Number of obs: 60, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.16891    0.03129  2.23939   5.398   0.0254 *
OriginTransplant -0.02845    0.01468 56.00000  -1.938   0.0577 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.235

eta_squared(TLE_T3.4.lme_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.06 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
TLE_T3.4.lme_SS<-lmer(log(TLE_cm.day+1)~Origin+(1|Genotype), data=Growth_T3.4[which(Growth_T3.4$Site=="SS"),])
summary(TLE_T3.4.lme_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin + (1 | Genotype)
   Data: Growth_T3.4[which(Growth_T3.4$Site == "SS"), ]

REML criterion at convergence: -185.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.44083 -0.56044  0.02132  0.52475  2.84423 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.001374 0.03707 
 Residual             0.001934 0.04398 
Number of obs: 60, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)       0.11023    0.02286  2.27122   4.822  0.03125 * 
OriginTransplant  0.03439    0.01135 56.00000   3.029  0.00371 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.248

eta_squared(TLE_T3.4.lme_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.14 | [0.03, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
TLE_T3.4.site.res<-data.frame(rbind(summary(TLE_T3.4.lme_KL)$coefficients[-1,],
                                     summary(TLE_T3.4.lme_SS)$coefficients[-1,]))
TLE_T3.4.site.res$Predictor<-c("KL Origin", "SS Origin")
TLE_T3.4.site.res$EtaSq<-c(eta_squared(TLE_T3.4.lme_KL)$Eta2, eta_squared(TLE_T3.4.lme_SS)$Eta2)
TLE_T3.4.site.res$Response<-rep("Growth", nrow(TLE_T3.4.site.res))
TLE_T3.4.site.res$Timepoint<-rep("TP3_4", nrow(TLE_T3.4.site.res))

##Combine results 
TLE_T3.4.res<-rbind(TLE_T3.4.res, TLE_T3.4.site.res)

##Summary statistics by Site and Origin
TP3.4_Growth.sum<-summarySE(Growth_T3.4, measurevar="TLE_cm.day", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Growth Rate across Treatments
TP3.4_Growth.plot<-ggplot(TP3.4_Growth.sum, aes(x=Site, y=TLE_cm.day, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=TLE_cm.day-se, ymax=TLE_cm.day+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y=expression(paste('Growth Rate (cm day'^-1*")")), colour="Origin")+
  ylim(0, 0.25)+
  annotate("text", x=1, y=0.24, label="-", size=sig.sz, fontface="bold")+
  annotate("text", x=2, y=0.2, label="**", size=sig.sz, fontface="bold")+
  geom_bracket(xmin=1, xmax=2, y.position=0.01, label.size=levels.sz, label="***", inherit.aes = FALSE); TP3.4_Growth.plot

Average Growth Rates by Genotype

SS

##Summary statistics by Genotype and Origin
TP3.4_Growth.sum.SS<-summarySE(subset(Growth_T3.4, Site=="SS"), measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin"), na.rm=TRUE)

TP3.4_Growth.sum.SS

##Calculate Percent Difference between Native and Transplant by Genotype
TP3.4_Growth.SS<-TP3.4_Growth.sum.SS[which(TP3.4_Growth.sum.SS$Origin=="Native"), c(1,4)]
names(TP3.4_Growth.SS)[2]<-paste(names(TP3.4_Growth.SS)[2], "N", sep="_")
TP3.4_Growth.SS$TLE_cm.day_T<-TP3.4_Growth.sum.SS$TLE_cm.day[which(TP3.4_Growth.sum.SS$Origin=="Transplant")]
TP3.4_Growth.SS$Perc.Inc<-((TP3.4_Growth.SS$TLE_cm.day_T-TP3.4_Growth.SS$TLE_cm.day_N)/TP3.4_Growth.SS$TLE_cm.day_N)*100

TP3.4_Growth.SS

##Mean and Standard Error
mean(TP3.4_Growth.SS$Perc.Inc)

[1] 34.22395

std.error(TP3.4_Growth.SS$Perc.Inc)

[1] 18.61694

KL

##Summary statistics by Genotype and Origin
TP3.4_Growth.sum.KL<-summarySE(subset(Growth_T3.4, Site=="KL"), measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin"), na.rm=TRUE)

TP3.4_Growth.sum.KL

##Calculate Percent Difference between Native and Transplant by Genotype
TP3.4_Growth.KL<-TP3.4_Growth.sum.KL[which(TP3.4_Growth.sum.KL$Origin=="Native"), c(1,4)]
names(TP3.4_Growth.KL)[2]<-paste(names(TP3.4_Growth.KL)[2], "N", sep="_")
TP3.4_Growth.KL$TLE_cm.day_T<-TP3.4_Growth.sum.KL$TLE_cm.day[which(TP3.4_Growth.sum.SS$Origin=="Transplant")]
TP3.4_Growth.KL$Perc.Inc<-((TP3.4_Growth.KL$TLE_cm.day_N-TP3.4_Growth.KL$TLE_cm.day_T)/TP3.4_Growth.KL$TLE_cm.day_T)*100

TP3.4_Growth.KL

##Mean and Standard Error
mean(TP3.4_Growth.KL$Perc.Inc)

[1] 21.37099

std.error(TP3.4_Growth.KL$Perc.Inc)

[1] 4.101183

Write Out Growth Data

##Growth Data dataframe
GrowthData<-Growth_T1.2

##Specify Timepoints
GrowthData<-GrowthData %>% dplyr::rename(T1.2_Ext_cm = Ext_cm, T1.2_TLE_cm.day = TLE_cm.day)

##Merge with Timepoints 3 to 4
GrowthData<-merge(GrowthData, Growth_T3.4, all=TRUE)

##Specify Timepoints
GrowthData<-GrowthData %>% dplyr::rename(T3.4_Ext_cm = Ext_cm, T3.4_TLE_cm.day = TLE_cm.day)

##Growth Data
write.csv(GrowthData, "Outputs/GrowthData.csv", row.names=FALSE)

Thermal Tolerance

Percent Retention

Subset Thermal Assay Data by Treatment

##Control
Therm_C<-subset(Thermal, Treat=="C")

##Heated
Therm_H<-subset(Thermal, Treat=="H")

Calculate Percent Retention

Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.

##Calculate averages for Control Treatment for each Site, Genotype, and Timepoint 
Therm_C<-na.omit(Therm_C)
names(Therm_C)

 [1] "ID"          "RandN"       "TimeP"       "Site"        "Genotype"    "Treat"      
 [7] "Treatment"   "Orig"        "Origin"      "Set"         "Site.Orig"   "SA_cm2"     
[13] "Chl_ug.cm2"  "Sym10.6_cm2" "Fv_Fm"      

Therm_C.a<-aggregate(Therm_C[,c(13:15)], list(Therm_C$Site, Therm_C$Genotype, Therm_C$TimeP, Therm_C$Origin), mean, na.action = na.omit)
names(Therm_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Origin")
names(Therm_C.a)[5:7]<-paste(names(Therm_C.a)[5:7], "C", sep="_")

##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, and Timepoint
Therm_H<-merge(Therm_H, Therm_C.a, all.x=TRUE )

##Calculate Proportion Retained for each Bleaching Metric relative to Control
Therm_H$Chl.prop<-round((Therm_H$Chl_ug.cm2/Therm_H$Chl_ug.cm2_C), 4)
Therm_H$Sym.prop<-round((Therm_H$Sym10.6_cm2/Therm_H$Sym10.6_cm2_C), 4)
Therm_H$PAM.prop<-round((Therm_H$Fv_Fm/Therm_H$Fv_Fm_C), 4)

##Set values >1 to 1
Therm_H$Chl.prop[which(Therm_H$Chl.prop>1)]<-1.0000
Therm_H$Sym.prop[which(Therm_H$Sym.prop>1)]<-1.0000
Therm_H$PAM.prop[which(Therm_H$PAM.prop>1)]<-1.0000

Initial Heat Assay

##Subset Initial Timepoint
Thermal_IN<-subset(Thermal, TimeP=="IN")
Thermal_IN$Site.Treat<-paste(Thermal_IN$Site, Thermal_IN$Treat, sep=".")
Thermal_IN$Site.Treat<-factor(Thermal_IN$Site.Treat, levels=c("KL.C", "KL.H", "SS.C", "SS.H"))


Therm_IN<-subset(Therm_H, TimeP=="IN")

Heated vs Control

Symbionts

#Check normality
hist(Thermal_IN$Sym10.6_cm2)

shapiro.test(Thermal_IN$Sym10.6_cm2)

    Shapiro-Wilk normality test

data:  Thermal_IN$Sym10.6_cm2
W = 0.90824, p-value = 0.001334

#Not normal

hist(log(Thermal_IN$Sym10.6_cm2+1))

shapiro.test(log(Thermal_IN$Sym10.6_cm2+1))

    Shapiro-Wilk normality test

data:  log(Thermal_IN$Sym10.6_cm2 + 1)
W = 0.95667, p-value = 0.07936

#Normal

##Model
#Function of Site and Treatment, with Genotype as a Random effect
Sym.lme_IN<-lmer(log(Sym10.6_cm2+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)

##Check residuals
Sym.lme_res_IN <- simulateResiduals(fittedModel = Sym.lme_IN, plot = F)
plot(Sym.lme_res_IN)

##Model results
summary(Sym.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym10.6_cm2 + 1) ~ Site * Treatment + (1 | Genotype)
   Data: Thermal_IN

REML criterion at convergence: -29.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.79325 -0.66013  0.08255  0.57806  2.12486 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0008222 0.02867 
 Residual             0.0230525 0.15183 
Number of obs: 47, groups:  Genotype, 3

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           0.38823    0.04685 14.06727   8.286 8.77e-07 ***
SiteSS                0.19508    0.06198 41.00308   3.147  0.00307 ** 
TreatmentHeat        -0.20370    0.06340 41.10985  -3.213  0.00256 ** 
SiteSS:TreatmentHeat -0.01724    0.08867 41.05799  -0.194  0.84677    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) SiteSS TrtmnH
SiteSS      -0.661              
TreatmentHt -0.647  0.489       
StSS:TrtmnH  0.462 -0.699 -0.715

eta_squared(Sym.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Site           |           0.30 | [0.12, 1.00]
Treatment      |           0.36 | [0.17, 1.00]
Site:Treatment |       9.20e-04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)

$emmeans
Treatment = Control:
 Site emmean     SE   df lower.CL upper.CL
 KL    0.388 0.0469 14.1   0.2878    0.489
 SS    0.583 0.0469 14.1   0.4829    0.684

Treatment = Heat:
 Site emmean     SE   df lower.CL upper.CL
 KL    0.185 0.0488 15.6   0.0808    0.288
 SS    0.362 0.0469 14.1   0.2619    0.463

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Treatment = Control:
 contrast estimate     SE   df t.ratio p.value
 KL - SS    -0.195 0.0620 41.0  -3.147  0.0031

Treatment = Heat:
 contrast estimate     SE   df t.ratio p.value
 KL - SS    -0.178 0.0635 41.1  -2.801  0.0077

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)

$emmeans
Site = KL:
 Treatment emmean     SE   df lower.CL upper.CL
 Control    0.388 0.0469 14.1   0.2878    0.489
 Heat       0.185 0.0488 15.6   0.0808    0.288

Site = SS:
 Treatment emmean     SE   df lower.CL upper.CL
 Control    0.583 0.0469 14.1   0.4829    0.684
 Heat       0.362 0.0469 14.1   0.2619    0.463

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast       estimate     SE   df t.ratio p.value
 Control - Heat    0.204 0.0635 41.1   3.209  0.0026

Site = SS:
 contrast       estimate     SE   df t.ratio p.value
 Control - Heat    0.221 0.0620 41.0   3.564  0.0009

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

##Save model results
Sym_IN.res<-data.frame(summary(Sym.lme_IN)$coefficients[-1,])

##Add pairwise
names(Sym_IN.res)<-names(data.frame(emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
Sym_IN.res<-rbind(Sym_IN.res,
data.frame(rbind(emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])

##Metadata
Sym_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
Sym_IN.res$Response<-rep("Symbionts", nrow(Sym_IN.res))
Sym_IN.res$EtaSq<-c(eta_squared(Sym.lme_IN)$Eta2, rep(NA, 4))

##Summary statistics by Site
IN_Sym.sum<-summarySE(Thermal_IN, measurevar="Sym10.6_cm2", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)

##Plot Average Symbionts across Treatments
IN_Sym.plot<-ggplot(IN_Sym.sum, aes(x=Site.Treat, y=Sym10.6_cm2, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Sym10.6_cm2-se, ymax=Sym10.6_cm2+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
  scale_shape_manual(values=c(1,16))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Treatment", y=expression(paste('Symbiont Density ('*10^6,'cells cm'^-2*")")), colour="Site")+
  ylim(0, 1); IN_Sym.plot

Chlorophyll

#Check normality
hist(Thermal_IN$Chl_ug.cm2)

shapiro.test(Thermal_IN$Chl_ug.cm2)

    Shapiro-Wilk normality test

data:  Thermal_IN$Chl_ug.cm2
W = 0.88059, p-value = 0.0002132

#Not normal

hist(log(Thermal_IN$Chl_ug.cm2+1))

shapiro.test(log(Thermal_IN$Chl_ug.cm2+1))

    Shapiro-Wilk normality test

data:  log(Thermal_IN$Chl_ug.cm2 + 1)
W = 0.92511, p-value = 0.005653

#Still not normal but less skewed

##Model
#Function of Site and Treatment, with Genotype as a Random effect
Chl.lme_IN<-lmer(log(Chl_ug.cm2+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)

##Check residuals
Chl.lme_res_IN <- simulateResiduals(fittedModel = Chl.lme_IN, plot = F)
plot(Chl.lme_res_IN)

##Model results
summary(Chl.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl_ug.cm2 + 1) ~ Site * Treatment + (1 | Genotype)
   Data: Thermal_IN

REML criterion at convergence: -69.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.50589 -0.51517  0.07352  0.50169  2.35158 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.006137 0.07834 
 Residual             0.007930 0.08905 
Number of obs: 46, groups:  Genotype, 3

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           0.45648    0.05202  2.95779   8.774  0.00329 ** 
SiteSS                0.27810    0.03635 39.99255   7.650 2.38e-09 ***
TreatmentHeat        -0.37295    0.03721 40.00578 -10.024 1.80e-12 ***
SiteSS:TreatmentHeat -0.16117    0.05265 40.01278  -3.061  0.00393 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) SiteSS TrtmnH
SiteSS      -0.349              
TreatmentHt -0.341  0.489       
StSS:TrtmnH  0.241 -0.691 -0.707

eta_squared(Chl.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Site           |           0.58 | [0.41, 1.00]
Treatment      |           0.88 | [0.82, 1.00]
Site:Treatment |           0.19 | [0.04, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)

$emmeans
Treatment = Control:
 Site emmean     SE   df lower.CL upper.CL
 KL   0.4565 0.0520 2.97   0.2899    0.623
 SS   0.7346 0.0520 2.97   0.5680    0.901

Treatment = Heat:
 Site emmean     SE   df lower.CL upper.CL
 KL   0.0835 0.0526 3.11  -0.0808    0.248
 SS   0.2005 0.0526 3.11   0.0361    0.365

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Treatment = Control:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.278 0.0364 40  -7.650  <.0001

Treatment = Heat:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.117 0.0381 40  -3.070  0.0038

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)

$emmeans
Site = KL:
 Treatment emmean     SE   df lower.CL upper.CL
 Control   0.4565 0.0520 2.97   0.2899    0.623
 Heat      0.0835 0.0526 3.11  -0.0808    0.248

Site = SS:
 Treatment emmean     SE   df lower.CL upper.CL
 Control   0.7346 0.0520 2.97   0.5680    0.901
 Heat      0.2005 0.0526 3.11   0.0361    0.365

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast       estimate     SE df t.ratio p.value
 Control - Heat    0.373 0.0372 40  10.022  <.0001

Site = SS:
 contrast       estimate     SE df t.ratio p.value
 Control - Heat    0.534 0.0372 40  14.353  <.0001

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

##Save model results
Chl_IN.res<-data.frame(summary(Chl.lme_IN)$coefficients[-1,])

##Add pairwise
names(Chl_IN.res)<-names(data.frame(emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
Chl_IN.res<-rbind(Chl_IN.res,
data.frame(rbind(emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])

##Metadata
Chl_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
Chl_IN.res$Response<-rep("Chlorophyll", nrow(Chl_IN.res))
Chl_IN.res$EtaSq<-c(eta_squared(Chl.lme_IN)$Eta2, rep(NA, 4))

##Summary statistics by Site
IN_Chl.sum<-summarySE(Thermal_IN, measurevar="Chl_ug.cm2", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)

##Plot Average Chlorophyll across Treatments
IN_Chl.plot<-ggplot(IN_Chl.sum, aes(x=Site.Treat, y=Chl_ug.cm2, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Chl_ug.cm2-se, ymax=Chl_ug.cm2+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
  scale_shape_manual(values=c(1,16))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Treatment", y=expression(paste('Total Chlorophyll (\u03BCg cm'^-2*")")), colour="Site")+
  ylim(0, 1.25); IN_Chl.plot

FvFm

#Check normality
hist(Thermal_IN$Fv_Fm)

shapiro.test(Thermal_IN$Fv_Fm)

    Shapiro-Wilk normality test

data:  Thermal_IN$Fv_Fm
W = 0.87006, p-value = 9.097e-05

#Not normal

hist(log(Thermal_IN$Fv_Fm+1))

shapiro.test(log(Thermal_IN$Fv_Fm+1))

    Shapiro-Wilk normality test

data:  log(Thermal_IN$Fv_Fm + 1)
W = 0.85365, p-value = 3.226e-05

#Still not normal but less skewed

##Model
#Function of Site and Treatment, with Genotype as a Random effect
PAM.lme_IN<-lmer(log(Fv_Fm+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)

boundary (singular) fit: see help('isSingular')

##Check residuals
PAM.lme_res_IN <- simulateResiduals(fittedModel = PAM.lme_IN, plot = F)
plot(PAM.lme_res_IN)

##Model results
summary(PAM.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Fv_Fm + 1) ~ Site * Treatment + (1 | Genotype)
   Data: Thermal_IN

REML criterion at convergence: -202.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7092 -0.4159  0.1391  0.6366  1.5720 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0000000 0.00000 
 Residual             0.0004174 0.02043 
Number of obs: 47, groups:  Genotype, 3

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           0.480202   0.005897 43.000000  81.425  < 2e-16 ***
SiteSS                0.013981   0.008340 43.000000   1.676    0.101    
TreatmentHeat        -0.044190   0.008528 43.000000  -5.182 5.57e-06 ***
SiteSS:TreatmentHeat  0.022162   0.011928 43.000000   1.858    0.070 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) SiteSS TrtmnH
SiteSS      -0.707              
TreatmentHt -0.692  0.489       
StSS:TrtmnH  0.494 -0.699 -0.715
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(PAM.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Site           |           0.29 | [0.12, 1.00]
Treatment      |           0.42 | [0.23, 1.00]
Site:Treatment |           0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)

$emmeans
Treatment = Control:
 Site emmean      SE   df lower.CL upper.CL
 KL    0.480 0.00590 21.6    0.468    0.492
 SS    0.494 0.00590 21.6    0.482    0.506

Treatment = Heat:
 Site emmean      SE   df lower.CL upper.CL
 KL    0.436 0.00618 23.5    0.423    0.449
 SS    0.472 0.00590 21.6    0.460    0.484

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Treatment = Control:
 contrast estimate      SE   df t.ratio p.value
 KL - SS   -0.0140 0.00834 41.0  -1.676  0.1013

Treatment = Heat:
 contrast estimate      SE   df t.ratio p.value
 KL - SS   -0.0361 0.00855 41.2  -4.229  0.0001

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)

$emmeans
Site = KL:
 Treatment emmean      SE   df lower.CL upper.CL
 Control    0.480 0.00590 21.6    0.468    0.492
 Heat       0.436 0.00618 23.5    0.423    0.449

Site = SS:
 Treatment emmean      SE   df lower.CL upper.CL
 Control    0.494 0.00590 21.6    0.482    0.506
 Heat       0.472 0.00590 21.6    0.460    0.484

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast       estimate      SE   df t.ratio p.value
 Control - Heat   0.0442 0.00855 41.2   5.171  <.0001

Site = SS:
 contrast       estimate      SE   df t.ratio p.value
 Control - Heat   0.0220 0.00834 41.0   2.641  0.0116

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

##Save model results
PAM_IN.res<-data.frame(summary(PAM.lme_IN)$coefficients[-1,])

##Add pairwise
names(PAM_IN.res)<-names(data.frame(emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
PAM_IN.res<-rbind(PAM_IN.res,
data.frame(rbind(emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])

##Metadata
PAM_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
PAM_IN.res$Response<-rep("FvFm", nrow(PAM_IN.res))
PAM_IN.res$EtaSq<-c(eta_squared(PAM.lme_IN)$Eta2, rep(NA, 4))

Not normal, but model residuals are OK.

##Summary statistics by Site
IN_PAM.sum<-summarySE(Thermal_IN, measurevar="Fv_Fm", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)

##Plot Average FvFm across Treatments
IN_PAM.plot<-ggplot(IN_PAM.sum, aes(x=Site.Treat, y=Fv_Fm, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Fv_Fm-se, ymax=Fv_Fm+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
  scale_shape_manual(values=c(1,16))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Treatment", y="Photochemical Efficiency (Fv/Fm)", colour="Site")+
  ylim(0.5, 0.655); IN_PAM.plot

Retention by Site

Symbionts

#Check normality
hist(Therm_IN$Sym.prop)

shapiro.test(Therm_IN$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_IN$Sym.prop
W = 0.93042, p-value = 0.1117

#Normal

##Model
#Function of Site, with Genotype as a Random effect
Tol_Sym.lme_IN<-lmer(Sym.prop~Site+(1|Genotype), data=Therm_IN)

##Check residuals
Tol_Sym.lme_res_IN <- simulateResiduals(fittedModel = Tol_Sym.lme_IN, plot = F)
plot(Tol_Sym.lme_res_IN)

##Model results
summary(Tol_Sym.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Site + (1 | Genotype)
   Data: Therm_IN

REML criterion at convergence: -2.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6802 -0.6578 -0.1521  0.6805  1.9852 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.07818  0.2796  
 Residual             0.03039  0.1743  
Number of obs: 23, groups:  Genotype, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)  0.46139    0.16983  2.21175   2.717   0.1014  
SiteSS       0.14316    0.07291 19.00328   1.964   0.0644 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
SiteSS -0.225

eta_squared(Tol_Sym.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Site      |           0.17 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_IN.res<-data.frame(summary(Tol_Sym.lme_IN)$coefficients)
names(Tol_Sym_IN.res)<-names(Sym_IN.res)[1:5]
Tol_Sym_IN.res$Predictor<-rep("Site", nrow(Tol_Sym_IN.res))
Tol_Sym_IN.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_IN.res))
Tol_Sym_IN.res$EtaSq<-c(eta_squared(Tol_Sym.lme_IN)$Eta2)

##Combine results
Sym_IN.res<-rbind(Sym_IN.res, Tol_Sym_IN.res[2,])

##Summary statistics by Site
IN_TolSym.sum<-summarySE(Therm_IN, measurevar="Sym.prop", groupvars=c("Site"), na.rm=TRUE)

##Plot Average Symbiont Retention across Treatments
IN_TolSym.plot<-ggplot(IN_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site", y="Symbiont Retention", colour="Site")+
  ylim(0, 1)+
  annotate("text", x=1.5, y=0.75, label="-", size=sig.sz, fontface="bold"); IN_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_IN$Chl.prop)

shapiro.test(Therm_IN$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_IN$Chl.prop
W = 0.98432, p-value = 0.9694

#Normal

##Model
#Function of Site, with Genotype as a Random effect
Tol_Chl.lme_IN<-lmer(Chl.prop~Site+(1|Genotype), data=Therm_IN)

##Check residuals
Tol_Chl.lme_res_IN <- simulateResiduals(fittedModel = Tol_Chl.lme_IN, plot = F)
plot(Tol_Chl.lme_res_IN)

##Model results
summary(Tol_Chl.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Chl.prop ~ Site + (1 | Genotype)
   Data: Therm_IN

REML criterion at convergence: -47.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3909 -0.5550 -0.1823  0.5284  1.9947 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.004576 0.06765 
 Residual             0.003436 0.05861 
Number of obs: 22, groups:  Genotype, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)  0.16552    0.04290  2.40867   3.858    0.045 *
SiteSS       0.04306    0.02514 18.06310   1.713    0.104  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
SiteSS -0.293

eta_squared(Tol_Chl.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Site      |           0.14 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_IN.res<-data.frame(summary(Tol_Chl.lme_IN)$coefficients)
names(Tol_Chl_IN.res)<-names(Chl_IN.res)[1:5]
Tol_Chl_IN.res$Predictor<-rep("Site", nrow(Tol_Chl_IN.res))
Tol_Chl_IN.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_IN.res))
Tol_Chl_IN.res$EtaSq<-c(eta_squared(Tol_Chl.lme_IN)$Eta2)

##Combine results
Chl_IN.res<-rbind(Chl_IN.res, Tol_Chl_IN.res[2,])

##Summary statistics by Site and Origin
IN_TolChl.sum<-summarySE(Therm_IN, measurevar="Chl.prop", groupvars=c("Site"), na.rm=TRUE)

##Plot Average Chlorophyll Retention across Treatments
IN_TolChl.plot<-ggplot(IN_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site", y="Chlorophyll Retention", colour="Site")+
  ylim(0, 0.35); IN_TolChl.plot

FvFm

#Check normality
hist(Therm_IN$PAM.prop)

shapiro.test(Therm_IN$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_IN$PAM.prop
W = 0.8289, p-value = 0.001154

#Not normal

hist(log(Therm_IN$PAM.prop+1))

shapiro.test(log(Therm_IN$PAM.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_IN$PAM.prop + 1)
W = 0.81061, p-value = 0.0005693

##Still not normal

##Model
#Function of Site, with Genotype as a Random effect
Tol_PAM.lme_IN<-lmer(PAM.prop~Site+(1|Genotype), data=Therm_IN)

##Check residuals
Tol_PAM.lme_res_IN <- simulateResiduals(fittedModel = Tol_PAM.lme_IN, plot = F)
plot(Tol_PAM.lme_res_IN)

##Model results
summary(Tol_PAM.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Site + (1 | Genotype)
   Data: Therm_IN

REML criterion at convergence: -50.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7860 -0.2498  0.1225  0.5362  1.4591 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002221 0.04713 
 Residual             0.003538 0.05948 
Number of obs: 23, groups:  Genotype, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.89138    0.03262  2.77551  27.325 0.000185 ***
SiteSS       0.05230    0.02487 19.00103   2.103 0.049053 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
SiteSS -0.399

eta_squared(Tol_PAM.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Site      |           0.19 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_IN.res<-data.frame(summary(Tol_PAM.lme_IN)$coefficients)
names(Tol_PAM_IN.res)<-names(PAM_IN.res)[1:5]
Tol_PAM_IN.res$Predictor<-rep("Site", nrow(Tol_PAM_IN.res))
Tol_PAM_IN.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_IN.res))
Tol_PAM_IN.res$EtaSq<-c(eta_squared(Tol_PAM.lme_IN)$Eta2)

##Combine results
PAM_IN.res<-rbind(PAM_IN.res, Tol_PAM_IN.res[2,])

Not normal, but model residuals are OK.

##Summary statistics by Site and Origin
IN_TolPAM.sum<-summarySE(Therm_IN, measurevar="PAM.prop", groupvars=c("Site"), na.rm=TRUE)

##Plot Average FvFm Retention across Treatments
IN_TolPAM.plot<-ggplot(IN_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site", y="Photochemical Efficiency Retention", colour="Site")+
  ylim(0.6, 1)+
  annotate("text", x=1.5, y=1, label="*", size=sig.sz, fontface="bold"); IN_TolPAM.plot

Tolerance after Transplant

TP1

##Subset Timepoint 1
Therm_TP1<-subset(Therm_H, TimeP=="TP1")

Symbionts

#Check normality
hist(Therm_TP1$Sym.prop)

shapiro.test(Therm_TP1$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP1$Sym.prop
W = 0.91699, p-value = 0.002332

#Not normal

hist(log(Therm_TP1$Sym.prop+1))

shapiro.test(log(Therm_TP1$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP1$Sym.prop + 1)
W = 0.92258, p-value = 0.003671

##Still not normal but less skewed

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP1<-lmer(log(Sym.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP1)

##Check residuals
Tol_Sym.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP1, plot = F)
plot(Tol_Sym.lme_res_TP1)

##Model results
summary(Tol_Sym.lme_TP1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP1

REML criterion at convergence: -33.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1715 -0.5714  0.0059  0.7146  1.9799 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.03376  0.1837  
 Residual             0.01886  0.1373  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)              0.27621    0.11324  2.42412   2.439   0.1128  
OriginTransplant         0.11668    0.05606 42.00000   2.081   0.0435 *
SiteSS                   0.06309    0.05606 42.00000   1.125   0.2668  
OriginTransplant:SiteSS -0.13250    0.07928 42.00000  -1.671   0.1021  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.248              
SiteSS      -0.248  0.500       
OrgnTrn:SSS  0.175 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP1)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |       1.51e-04 | [0.00, 1.00]
Origin:Site |           0.06 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP1.res<-data.frame(summary(Tol_Sym.lme_TP1)$coefficients[-1,])
Tol_Sym_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP1.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP1)$Eta2)

Not normal, but model residuals are OK.

Effect size of Origin for each Site

##KL
Tol_Sym.lme_TP1_KL<-lmer(log(Sym.prop+1)~Origin+(1|Genotype), data=Therm_TP1[which(Therm_TP1$Site=="KL"),])
summary(Tol_Sym.lme_TP1_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP1[which(Therm_TP1$Site == "KL"), ]

REML criterion at convergence: -21.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.22211 -0.64466  0.06469  0.75749  1.28798 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.04259  0.2064  
 Residual             0.01291  0.1136  
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.27621    0.12358  2.14836   2.235   0.1461  
OriginTransplant  0.11668    0.04638 20.00000   2.516   0.0205 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.188

eta_squared(Tol_Sym.lme_TP1_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.24 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Sym.lme_TP1_SS<-lmer(log(Sym.prop+1)~Origin+(1|Genotype), data=Therm_TP1[which(Therm_TP1$Site=="SS"),])
summary(Tol_Sym.lme_TP1_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP1[which(Therm_TP1$Site == "SS"), ]

REML criterion at convergence: -8.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4670 -0.5662  0.1529  0.6190  1.4817 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.02355  0.1535  
 Residual             0.02580  0.1606  
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.33930    0.10000  2.50724   3.393   0.0557 .
OriginTransplant -0.01583    0.06558 20.00000  -0.241   0.8118  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.328

eta_squared(Tol_Sym.lme_TP1_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |       2.90e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Sym_TP1.site.res<-data.frame(rbind(summary(Tol_Sym.lme_TP1_KL)$coefficients[-1,],
                                     summary(Tol_Sym.lme_TP1_SS)$coefficients[-1,]))
Tol_Sym_TP1.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Sym_TP1.site.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP1_KL)$Eta2, eta_squared(Tol_Sym.lme_TP1_SS)$Eta2)

##Combine results 
Tol_Sym_TP1.res<-rbind(Tol_Sym_TP1.res, Tol_Sym_TP1.site.res)
Tol_Sym_TP1.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP1.res))

##Summary statistics by Site and Origin
TP1_TolSym.sum<-summarySE(Therm_TP1, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Symbiont Retention across Treatments
TP1_TolSym.plot<-ggplot(TP1_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1)+
  annotate("text", x=1, y=0.75, label="*", size=sig.sz, fontface="bold"); TP1_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP1$Chl.prop)

shapiro.test(Therm_TP1$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP1$Chl.prop
W = 0.94139, p-value = 0.01836

#Not normal

hist(log(Therm_TP1$Chl.prop+1))

shapiro.test(log(Therm_TP1$Chl.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP1$Chl.prop + 1)
W = 0.94765, p-value = 0.03227

##Still not normal but less skewed

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP1<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP1)

##Check residuals
Tol_Chl.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP1, plot = F)
plot(Tol_Chl.lme_res_TP1)

##Model results
summary(Tol_Chl.lme_TP1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP1

REML criterion at convergence: -112.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1439 -0.5458 -0.0382  0.6267  1.8760 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.008382 0.09155 
 Residual             0.003065 0.05537 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)
(Intercept)              0.13931    0.05522  2.27664   2.523    0.113
OriginTransplant         0.01521    0.02260 42.00000   0.673    0.505
SiteSS                   0.02628    0.02260 42.00000   1.163    0.251
OriginTransplant:SiteSS -0.02689    0.03197 42.00000  -0.841    0.405

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.205              
SiteSS      -0.205  0.500       
OrgnTrn:SSS  0.145 -0.707 -0.707

eta_squared(Tol_Chl.lme_TP1)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       2.90e-04 | [0.00, 1.00]
Site        |           0.02 | [0.00, 1.00]
Origin:Site |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP1.res<-data.frame(summary(Tol_Chl.lme_TP1)$coefficients[-1,])
Tol_Chl_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP1.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP1)$Eta2)

Effect size of Origin for each Site

##KL
Tol_Chl.lme_TP1_KL<-lmer(log(Chl.prop+1)~Origin+(1|Genotype), data=Therm_TP1[which(Therm_TP1$Site=="KL"),])
summary(Tol_Chl.lme_TP1_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP1[which(Therm_TP1$Site == "KL"), ]

REML criterion at convergence: -55.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.04878 -0.58286  0.01887  0.57011  1.44796 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.007318 0.08555 
 Residual             0.002846 0.05334 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)
(Intercept)       0.13931    0.05173  2.18924   2.693    0.104
OriginTransplant  0.01521    0.02178 20.00000   0.698    0.493

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.210

eta_squared(Tol_Chl.lme_TP1_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Chl.lme_TP1_SS<-lmer(log(Chl.prop+1)~Origin+(1|Genotype), data=Therm_TP1[which(Therm_TP1$Site=="SS"),])
summary(Tol_Chl.lme_TP1_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP1[which(Therm_TP1$Site == "SS"), ]

REML criterion at convergence: -50.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.43149 -0.56754 -0.06902  0.58853  1.85344 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.009135 0.09558 
 Residual             0.003511 0.05926 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.16560    0.05777  2.18712   2.866   0.0931 .
OriginTransplant -0.01168    0.02419 20.00000  -0.483   0.6344  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.209

eta_squared(Tol_Chl.lme_TP1_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.01 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Chl_TP1.site.res<-data.frame(rbind(summary(Tol_Chl.lme_TP1_KL)$coefficients[-1,],
                                     summary(Tol_Chl.lme_TP1_SS)$coefficients[-1,]))
Tol_Chl_TP1.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Chl_TP1.site.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP1_KL)$Eta2, eta_squared(Tol_Chl.lme_TP1_SS)$Eta2)

##Combine results 
Tol_Chl_TP1.res<-rbind(Tol_Chl_TP1.res, Tol_Chl_TP1.site.res)
Tol_Chl_TP1.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP1.res))

##Summary statistics by Site and Origin
TP1_TolChl.sum<-summarySE(Therm_TP1, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Chlorophyll Retention across Treatments
TP1_TolChl.plot<-ggplot(TP1_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 0.35); TP1_TolChl.plot

FvFm

#Check normality
hist(Therm_TP1$PAM.prop)

shapiro.test(Therm_TP1$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP1$PAM.prop
W = 0.84073, p-value = 1.245e-05

#Not normal

hist(log(Therm_TP1$PAM.prop+1))

shapiro.test(log(Therm_TP1$PAM.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP1$PAM.prop + 1)
W = 0.82049, p-value = 3.908e-06

##Still not normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP1<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP1)

##Check residuals
Tol_PAM.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP1, plot = F)
plot(Tol_PAM.lme_res_TP1)

##Model results
summary(Tol_PAM.lme_TP1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP1

REML criterion at convergence: -84.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3058 -0.3498  0.1733  0.6717  1.7818 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002735 0.05230 
 Residual             0.006214 0.07883 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.86181    0.03781  3.74555  22.794 3.72e-05 ***
OriginTransplant         0.06813    0.03218 42.00000   2.117   0.0402 *  
SiteSS                  -0.01277    0.03218 42.00000  -0.397   0.6936    
OriginTransplant:SiteSS -0.07942    0.04551 42.00000  -1.745   0.0883 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.426              
SiteSS      -0.426  0.500       
OrgnTrn:SSS  0.301 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP1)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |           0.11 | [0.01, 1.00]
Origin:Site |           0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP1.res<-data.frame(summary(Tol_PAM.lme_TP1)$coefficients[-1,])
Tol_PAM_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP1.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP1)$Eta2)

Not normal, but model residuals are OK.

Effect size of Origin for each Site

##KL
Tol_PAM.lme_TP1_KL<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP1[which(Therm_TP1$Site=="KL"),])
summary(Tol_PAM.lme_TP1_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP1[which(Therm_TP1$Site == "KL"), ]

REML criterion at convergence: -57.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.21445 -0.36299  0.09321  0.61667  1.51097 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0009092 0.03015 
 Residual             0.0031116 0.05578 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.86181    0.02371  3.34791  36.341 1.82e-05 ***
OriginTransplant  0.06813    0.02277 20.00000   2.992  0.00721 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.480

eta_squared(Tol_PAM.lme_TP1_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.31 | [0.06, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_PAM.lme_TP1_SS<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP1[which(Therm_TP1$Site=="SS"),])
summary(Tol_PAM.lme_TP1_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP1[which(Therm_TP1$Site == "SS"), ]

REML criterion at convergence: -33.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.69510 -0.36658  0.09367  0.51233  1.74473 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.005306 0.07284 
 Residual             0.008775 0.09367 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.84904    0.05000  2.73586  16.981 0.000746 ***
OriginTransplant -0.01128    0.03824 20.00000  -0.295 0.770994    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.382

eta_squared(Tol_PAM.lme_TP1_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |       4.33e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_PAM_TP1.site.res<-data.frame(rbind(summary(Tol_PAM.lme_TP1_KL)$coefficients[-1,],
                                     summary(Tol_PAM.lme_TP1_SS)$coefficients[-1,]))
Tol_PAM_TP1.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_PAM_TP1.site.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP1_KL)$Eta2, eta_squared(Tol_PAM.lme_TP1_SS)$Eta2)

##Combine results 
Tol_PAM_TP1.res<-rbind(Tol_PAM_TP1.res, Tol_PAM_TP1.site.res)
Tol_PAM_TP1.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP1.res))

##Summary statistics by Site and Origin
TP1_TolPAM.sum<-summarySE(Therm_TP1, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average FvFm Retention across Treatments
TP1_TolPAM.plot<-ggplot(TP1_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.6, 1)+
  annotate("text", x=1, y=1, label="**", size=sig.sz, fontface="bold"); TP1_TolPAM.plot

Save Results

##Combine Results
Tol_T1.res<-rbind(Tol_PAM_TP1.res, Tol_Sym_TP1.res, Tol_Chl_TP1.res)

##Add Timepoint
Tol_T1.res$TimeP<-rep("TP1", nrow(Tol_T1.res))

TP2

##Subset Timepoint 2
Therm_TP2<-subset(Therm_H, TimeP=="TP2")

Symbionts

#Check normality
hist(Therm_TP2$Sym.prop)

shapiro.test(Therm_TP2$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP2$Sym.prop
W = 0.89527, p-value = 0.0004439

#Not normal

hist(log(Therm_TP2$Sym.prop+1))

shapiro.test(log(Therm_TP2$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP2$Sym.prop + 1)
W = 0.88978, p-value = 0.000299

##Still not normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP2<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP2)

##Check residuals
Tol_Sym.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP2, plot = F)
plot(Tol_Sym.lme_res_TP2)

##Model results
summary(Tol_Sym.lme_TP2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP2

REML criterion at convergence: -3.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.20381 -0.54227  0.01202  0.47665  2.09092 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.04681  0.2164  
 Residual             0.03767  0.1941  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)              0.805833   0.136905  2.613359   5.886   0.0142 *
OriginTransplant        -0.107900   0.079237 42.000000  -1.362   0.1805  
SiteSS                  -0.178425   0.079237 42.000000  -2.252   0.0296 *
OriginTransplant:SiteSS  0.005183   0.112058 42.000000   0.046   0.9633  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.289              
SiteSS      -0.289  0.500       
OrgnTrn:SSS  0.205 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP2)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.08 | [0.00, 1.00]
Site        |           0.19 | [0.04, 1.00]
Origin:Site |       5.09e-05 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP2.res<-data.frame(summary(Tol_Sym.lme_TP2)$coefficients[-1,])
Tol_Sym_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP2.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP2)$Eta2)

Not normal, but model residuals are OK.

Effect size of Origin for each Site

##KL
Tol_Sym.lme_TP2_KL<-lmer(Sym.prop~Origin+(1|Genotype), data=Therm_TP2[which(Therm_TP2$Site=="KL"),])
summary(Tol_Sym.lme_TP2_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP2[which(Therm_TP2$Site == "KL"), ]

REML criterion at convergence: -1.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0500 -0.8276  0.2039  0.4804  2.0614 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.04096  0.2024  
 Residual             0.03487  0.1867  
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.80583    0.12868  2.40096   6.262   0.0154 *
OriginTransplant -0.10790    0.07623 20.00000  -1.415   0.1723  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.296

eta_squared(Tol_Sym.lme_TP2_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.09 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Sym.lme_TP2_SS<-lmer(Sym.prop~Origin+(1|Genotype), data=Therm_TP2[which(Therm_TP2$Site=="SS"),])
summary(Tol_Sym.lme_TP2_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP2[which(Therm_TP2$Site == "SS"), ]

REML criterion at convergence: 3.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.98373 -0.46709 -0.03333  0.50742  1.82777 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.04792  0.2189  
 Residual             0.04392  0.2096  
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.62741    0.14011  2.42980   4.478    0.032 *
OriginTransplant -0.10272    0.08556 20.00000  -1.201    0.244  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.305

eta_squared(Tol_Sym.lme_TP2_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Sym_TP2.site.res<-data.frame(rbind(summary(Tol_Sym.lme_TP2_KL)$coefficients[-1,],
                                     summary(Tol_Sym.lme_TP2_SS)$coefficients[-1,]))
Tol_Sym_TP2.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Sym_TP2.site.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP2_KL)$Eta2, eta_squared(Tol_Sym.lme_TP2_SS)$Eta2)

##Combine results 
Tol_Sym_TP2.res<-rbind(Tol_Sym_TP2.res, Tol_Sym_TP2.site.res)
Tol_Sym_TP2.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP2.res))

##Summary statistics by Site and Origin
TP2_TolSym.sum<-summarySE(Therm_TP2, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP2_TolSym.plot<-ggplot(TP2_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1); TP2_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP2$Chl.prop)

shapiro.test(Therm_TP2$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP2$Chl.prop
W = 0.97677, p-value = 0.4523

#Normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP2<-lmer(Chl.prop~Origin*Site+(1|Genotype), data=Therm_TP2)

##Check residuals
Tol_Chl.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP2, plot = F)
plot(Tol_Chl.lme_res_TP2)

##Model results
summary(Tol_Chl.lme_TP2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Chl.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP2

REML criterion at convergence: -81

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0024 -0.7622 -0.1174  0.6824  2.9666 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.006146 0.07840 
 Residual             0.006538 0.08086 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)              0.31077    0.05093  2.81331   6.102   0.0106 *
OriginTransplant        -0.03844    0.03301 42.00000  -1.165   0.2508  
SiteSS                  -0.06514    0.03301 42.00000  -1.973   0.0550 .
OriginTransplant:SiteSS  0.06166    0.04668 42.00000   1.321   0.1937  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.324              
SiteSS      -0.324  0.500       
OrgnTrn:SSS  0.229 -0.707 -0.707

eta_squared(Tol_Chl.lme_TP2)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       2.53e-03 | [0.00, 1.00]
Site        |           0.05 | [0.00, 1.00]
Origin:Site |           0.04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP2.res<-data.frame(summary(Tol_Chl.lme_TP2)$coefficients[-1,])
Tol_Chl_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP2.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP2)$Eta2)

Effect size of Origin for each Site

##KL
Tol_Chl.lme_TP2_KL<-lmer(Chl.prop~Origin+(1|Genotype), data=Therm_TP2[which(Therm_TP2$Site=="KL"),])
summary(Tol_Chl.lme_TP2_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Chl.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP2[which(Therm_TP2$Site == "KL"), ]

REML criterion at convergence: -39.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.9725 -0.7459 -0.3411  0.6270  2.6035 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.006040 0.07772 
 Residual             0.006269 0.07918 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.31077    0.05036  2.48246   6.171   0.0146 *
OriginTransplant -0.03844    0.03232 20.00001  -1.189   0.2483  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.321

eta_squared(Tol_Chl.lme_TP2_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Chl.lme_TP2_SS<-lmer(Chl.prop~Origin+(1|Genotype), data=Therm_TP2[which(Therm_TP2$Site=="SS"),])
summary(Tol_Chl.lme_TP2_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Chl.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP2[which(Therm_TP2$Site == "SS"), ]

REML criterion at convergence: -46

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2119 -0.6418  0.2701  0.7660  1.5135 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.009505 0.09749 
 Residual             0.004441 0.06664 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)       0.24563    0.05948  2.22611   4.129   0.0448 *
OriginTransplant  0.02322    0.02720 20.00000   0.853   0.4035  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.229

eta_squared(Tol_Chl.lme_TP2_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Chl_TP2.site.res<-data.frame(rbind(summary(Tol_Chl.lme_TP2_KL)$coefficients[-1,],
                                     summary(Tol_Chl.lme_TP2_SS)$coefficients[-1,]))
Tol_Chl_TP2.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Chl_TP2.site.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP2_KL)$Eta2, eta_squared(Tol_Chl.lme_TP2_SS)$Eta2)

##Combine results 
Tol_Chl_TP2.res<-rbind(Tol_Chl_TP2.res, Tol_Chl_TP2.site.res)
Tol_Chl_TP2.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP2.res))

##Summary statistics by Site and Origin
TP2_TolChl.sum<-summarySE(Therm_TP2, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP2_TolChl.plot<-ggplot(TP2_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 0.35); TP2_TolChl.plot

FvFm

#Check normality
hist(Therm_TP2$PAM.prop)

shapiro.test(Therm_TP2$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP2$PAM.prop
W = 0.73527, p-value = 5.896e-08

#Not normal

hist(log(Therm_TP2$PAM.prop+1))

shapiro.test(log(Therm_TP2$PAM.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP2$PAM.prop + 1)
W = 0.70205, p-value = 1.443e-08

##Still not normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP2<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP2)

##Check residuals
Tol_PAM.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP2, plot = F)
plot(Tol_PAM.lme_res_TP2)

##Model results
summary(Tol_PAM.lme_TP2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP2

REML criterion at convergence: -102.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8326 -0.3753 -0.0042  0.5350  1.7496 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.001675 0.04093 
 Residual             0.004192 0.06475 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.85622    0.03013  3.92146  28.419 1.09e-05 ***
OriginTransplant         0.07526    0.02643 42.00000   2.847  0.00680 ** 
SiteSS                   0.06487    0.02643 42.00000   2.454  0.01834 *  
OriginTransplant:SiteSS -0.10206    0.03738 42.00000  -2.730  0.00921 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.439              
SiteSS      -0.439  0.500       
OrgnTrn:SSS  0.310 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP2)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |           0.01 | [0.00, 1.00]
Origin:Site |           0.15 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP2.res<-data.frame(summary(Tol_PAM.lme_TP2)$coefficients[-1,])
Tol_PAM_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP2.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP2)$Eta2)

Not normal, but model residuals are OK.

Effect size of Origin for each Site

##KL
Tol_PAM.lme_TP2_KL<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP2[which(Therm_TP2$Site=="KL"),])
summary(Tol_PAM.lme_TP2_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP2[which(Therm_TP2$Site == "KL"), ]

REML criterion at convergence: -40.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9208 -0.4020  0.2520  0.5255  1.0060 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.003195 0.05653 
 Residual             0.006508 0.08067 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.85622    0.04009  2.88245  21.356  0.00029 ***
OriginTransplant  0.07526    0.03293 19.99994   2.285  0.03336 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.411

eta_squared(Tol_PAM.lme_TP2_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.21 | [0.01, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_PAM.lme_TP2_SS<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP2[which(Therm_TP2$Site=="SS"),])
summary(Tol_PAM.lme_TP2_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP2[which(Therm_TP2$Site == "SS"), ]

REML criterion at convergence: -71.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.04910 -0.74292 -0.04708  0.53594  2.26126 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0004937 0.02222 
 Residual             0.0016296 0.04037 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.92109    0.01733  3.31062  53.146 5.72e-06 ***
OriginTransplant -0.02680    0.01648 20.00000  -1.626     0.12    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.475

eta_squared(Tol_PAM.lme_TP2_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.12 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_PAM_TP2.site.res<-data.frame(rbind(summary(Tol_PAM.lme_TP2_KL)$coefficients[-1,],
                                     summary(Tol_PAM.lme_TP2_SS)$coefficients[-1,]))
Tol_PAM_TP2.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_PAM_TP2.site.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP2_KL)$Eta2, eta_squared(Tol_PAM.lme_TP2_SS)$Eta2)

##Combine results 
Tol_PAM_TP2.res<-rbind(Tol_PAM_TP2.res, Tol_PAM_TP2.site.res)
Tol_PAM_TP2.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP2.res))

##Summary statistics by Site and Origin
TP2_TolPAM.sum<-summarySE(Therm_TP2, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP2_TolPAM.plot<-ggplot(TP2_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.6, 1)+
  annotate("text", x=1, y=1, label="*", size=sig.sz, fontface="bold"); TP2_TolPAM.plot

Save Results

##Combine Results
Tol_T2.res<-rbind(Tol_PAM_TP2.res, Tol_Sym_TP2.res, Tol_Chl_TP2.res)

##Add Timepoint
Tol_T2.res$TimeP<-rep("TP2", nrow(Tol_T2.res))

TP3

##Subset Timepoint 3
Therm_TP3<-subset(Therm_H, TimeP=="TP3")

Symbionts

#Check normality
hist(Therm_TP3$Sym.prop)

shapiro.test(Therm_TP3$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP3$Sym.prop
W = 0.78653, p-value = 6.531e-07

#Not normal

hist(log(Therm_TP3$Sym.prop+1))

shapiro.test(log(Therm_TP3$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP3$Sym.prop + 1)
W = 0.78571, p-value = 6.267e-07

##Still not normal but less skewed

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP3<-lmer(log(Sym.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP3)

##Check residuals
Tol_Sym.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP3, plot = F)
plot(Tol_Sym.lme_res_TP3)

##Model results
summary(Tol_Sym.lme_TP3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP3

REML criterion at convergence: -108.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1559 -0.5389  0.3778  0.5963  1.6051 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 8.41e-05 0.00917 
 Residual             3.90e-03 0.06245 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.59780    0.01879 16.87446  31.815   <2e-16 ***
OriginTransplant         0.03978    0.02550 42.00000   1.560   0.1262    
SiteSS                   0.05892    0.02550 42.00000   2.311   0.0258 *  
OriginTransplant:SiteSS -0.03130    0.03606 42.00000  -0.868   0.3903    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.678              
SiteSS      -0.678  0.500       
OrgnTrn:SSS  0.480 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP3)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |           0.12 | [0.01, 1.00]
Origin:Site |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP3.res<-data.frame(summary(Tol_Sym.lme_TP3)$coefficients[-1,])
Tol_Sym_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP3.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP3)$Eta2)

Not normal, but model residuals are OK.

Effect size of Origin for each Site

##KL
Tol_Sym.lme_TP3_KL<-lmer(log(Sym.prop+1)~Origin+(1|Genotype), data=Therm_TP3[which(Therm_TP3$Site=="KL"),])
summary(Tol_Sym.lme_TP3_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP3[which(Therm_TP3$Site == "KL"), ]

REML criterion at convergence: -48.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.66523 -0.53461  0.01898  0.71797  1.50619 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0003833 0.01958 
 Residual             0.0048325 0.06952 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.59780    0.02303  5.00842  25.955 1.56e-06 ***
OriginTransplant  0.03978    0.02838 20.00000   1.402    0.176    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.616

eta_squared(Tol_Sym.lme_TP3_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.09 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Sym.lme_TP3_SS<-lmer(log(Sym.prop+1)~Origin+(1|Genotype), data=Therm_TP3[which(Therm_TP3$Site=="SS"),])

boundary (singular) fit: see help('isSingular')

summary(Tol_Sym.lme_TP3_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP3[which(Therm_TP3$Site == "SS"), ]

REML criterion at convergence: -61.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6314 -0.4640  0.5270  0.5670  0.6869 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.000000 0.00000 
 Residual             0.002812 0.05302 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.656722   0.015307 22.000000  42.904   <2e-16 ***
OriginTransplant  0.008483   0.021647 22.000000   0.392    0.699    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Sym.lme_TP3_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |       6.93e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Sym_TP3.site.res<-data.frame(rbind(summary(Tol_Sym.lme_TP3_KL)$coefficients[-1,],
                                     summary(Tol_Sym.lme_TP3_SS)$coefficients[-1,]))
Tol_Sym_TP3.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Sym_TP3.site.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP3_KL)$Eta2, eta_squared(Tol_Sym.lme_TP3_SS)$Eta2)

##Combine results 
Tol_Sym_TP3.res<-rbind(Tol_Sym_TP3.res, Tol_Sym_TP3.site.res)
Tol_Sym_TP3.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP3.res))

##Summary statistics by Site and Origin
TP3_TolSym.sum<-summarySE(Therm_TP3, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP3_TolSym.plot<-ggplot(TP3_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1); TP3_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP3$Chl.prop)

shapiro.test(Therm_TP3$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP3$Chl.prop
W = 0.89534, p-value = 0.0004461

#Not normal

hist(log(Therm_TP3$Chl.prop+1))

shapiro.test(log(Therm_TP3$Chl.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP3$Chl.prop + 1)
W = 0.9211, p-value = 0.003251

##Still not normal but less skewed

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP3<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP3)

boundary (singular) fit: see help('isSingular')

##Check residuals
Tol_Chl.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP3, plot = F)
plot(Tol_Chl.lme_res_TP3)

##Model results
summary(Tol_Chl.lme_TP3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP3

REML criterion at convergence: -47.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7095 -0.6080 -0.1405  0.4999  1.7650 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.00000  0.0000  
 Residual             0.01594  0.1262  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.43851    0.03644 44.00000  12.032 1.65e-15 ***
OriginTransplant         0.03181    0.05154 44.00000   0.617    0.540    
SiteSS                   0.03852    0.05154 44.00000   0.747    0.459    
OriginTransplant:SiteSS -0.02511    0.07289 44.00000  -0.344    0.732    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.707              
SiteSS      -0.707  0.500       
OrgnTrn:SSS  0.500 -0.707 -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP3)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       6.31e-03 | [0.00, 1.00]
Site        |           0.01 | [0.00, 1.00]
Origin:Site |       2.69e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP3.res<-data.frame(summary(Tol_Chl.lme_TP3)$coefficients[-1,])
Tol_Chl_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP3.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP3)$Eta2)

Not normal, but model residuals are OK.

Effect size of Origin for each Site

##KL
Tol_Chl.lme_TP3_KL<-lmer(log(Chl.prop+1)~Origin+(1|Genotype), data=Therm_TP3[which(Therm_TP3$Site=="KL"),])
summary(Tol_Chl.lme_TP3_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP3[which(Therm_TP3$Site == "KL"), ]

REML criterion at convergence: -35

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1966 -0.5212 -0.1283  0.4729  1.9007 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002703 0.05199 
 Residual             0.008473 0.09205 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)       0.43851    0.04009  3.25876  10.938  0.00112 **
OriginTransplant  0.03181    0.03758 20.00000   0.847  0.40725   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.469

eta_squared(Tol_Chl.lme_TP3_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Chl.lme_TP3_SS<-lmer(log(Chl.prop+1)~Origin+(1|Genotype), data=Therm_TP3[which(Therm_TP3$Site=="SS"),])

boundary (singular) fit: see help('isSingular')

summary(Tol_Chl.lme_TP3_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP3[which(Therm_TP3$Site == "SS"), ]

REML criterion at convergence: -17.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4740 -0.7462 -0.2782  0.6974  1.4761 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 Genotype (Intercept) 3.210e-22 1.792e-11
 Residual             2.144e-02 1.464e-01
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.477025   0.042265 22.000000  11.286 1.28e-10 ***
OriginTransplant  0.006709   0.059772 22.000000   0.112    0.912    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP3_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |       5.72e-04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Chl_TP3.site.res<-data.frame(rbind(summary(Tol_Chl.lme_TP3_KL)$coefficients[-1,],
                                     summary(Tol_Chl.lme_TP3_SS)$coefficients[-1,]))
Tol_Chl_TP3.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Chl_TP3.site.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP3_KL)$Eta2, eta_squared(Tol_Chl.lme_TP3_SS)$Eta2)

##Combine results 
Tol_Chl_TP3.res<-rbind(Tol_Chl_TP3.res, Tol_Chl_TP3.site.res)
Tol_Chl_TP3.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP3.res))

##Summary statistics by Site and Origin
TP3_TolChl.sum<-summarySE(Therm_TP3, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP3_TolChl.plot<-ggplot(TP3_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 1); TP3_TolChl.plot

FvFm

#Check normality
hist(Therm_TP3$PAM.prop)

shapiro.test(Therm_TP3$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP3$PAM.prop
W = 0.96358, p-value = 0.1408

#Normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP3<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP3)

##Check residuals
Tol_PAM.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP3, plot = F)
plot(Tol_PAM.lme_res_TP3)

##Model results
summary(Tol_PAM.lme_TP3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP3

REML criterion at convergence: -172.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.29399 -0.71463  0.01706  0.83761  1.77082 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 3.177e-05 0.005637
 Residual             9.061e-04 0.030101
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)              0.954325   0.009279 14.520200 102.850   <2e-16 ***
OriginTransplant         0.009025   0.012289 42.000000   0.734    0.467    
SiteSS                   0.003975   0.012289 42.000000   0.323    0.748    
OriginTransplant:SiteSS -0.019133   0.017379 42.000000  -1.101    0.277    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.662              
SiteSS      -0.662  0.500       
OrgnTrn:SSS  0.468 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP3)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       9.25e-05 | [0.00, 1.00]
Site        |       9.76e-03 | [0.00, 1.00]
Origin:Site |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP3.res<-data.frame(summary(Tol_PAM.lme_TP3)$coefficients[-1,])
Tol_PAM_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP3.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP3)$Eta2)

Effect size of Origin for each Site

##KL
Tol_PAM.lme_TP3_KL<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP3[which(Therm_TP3$Site=="KL"),])
summary(Tol_PAM.lme_TP3_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP3[which(Therm_TP3$Site == "KL"), ]

REML criterion at convergence: -87.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4364 -0.5727  0.1072  0.8113  1.2811 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 3.996e-05 0.006322
 Residual             8.452e-04 0.029072
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.954325   0.009152  5.656770  104.28 1.65e-10 ***
OriginTransplant  0.009025   0.011869 20.001145    0.76    0.456    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.648

eta_squared(Tol_PAM.lme_TP3_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_PAM.lme_TP3_SS<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP3[which(Therm_TP3$Site=="SS"),])

boundary (singular) fit: see help('isSingular')

summary(Tol_PAM.lme_TP3_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP3[which(Therm_TP3$Site == "SS"), ]

REML criterion at convergence: -84.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6542 -0.7178 -0.1194  0.6304  1.6515 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0000000 0.00000 
 Residual             0.0009841 0.03137 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.958300   0.009056 22.000000 105.820   <2e-16 ***
OriginTransplant -0.010108   0.012807 22.000000  -0.789    0.438    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_PAM.lme_TP3_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_PAM_TP3.site.res<-data.frame(rbind(summary(Tol_PAM.lme_TP3_KL)$coefficients[-1,],
                                     summary(Tol_PAM.lme_TP3_SS)$coefficients[-1,]))
Tol_PAM_TP3.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_PAM_TP3.site.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP3_KL)$Eta2, eta_squared(Tol_PAM.lme_TP3_SS)$Eta2)

##Combine results 
Tol_PAM_TP3.res<-rbind(Tol_PAM_TP3.res, Tol_PAM_TP3.site.res)
Tol_PAM_TP3.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP3.res))

##Summary statistics by Site and Origin
TP3_TolPAM.sum<-summarySE(Therm_TP3, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP3_TolPAM.plot<-ggplot(TP3_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.7, 1); TP3_TolPAM.plot

Save Results

##Combine Results
Tol_T3.res<-rbind(Tol_PAM_TP3.res, Tol_Sym_TP3.res, Tol_Chl_TP3.res)

##Add Timepoint
Tol_T3.res$TimeP<-rep("TP3", nrow(Tol_T3.res))

TP4

##Subset Timepoint 4
Therm_TP4<-subset(Therm_H, TimeP=="TP4")

Symbionts

#Check normality
hist(Therm_TP4$Sym.prop)

shapiro.test(Therm_TP4$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP4$Sym.prop
W = 0.93453, p-value = 0.01005

#Not normal

hist(log(Therm_TP4$Sym.prop+1))

shapiro.test(log(Therm_TP4$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP4$Sym.prop + 1)
W = 0.95458, p-value = 0.061

#Normal

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP4<-lmer(log(Sym.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP4)

##Check residuals
Tol_Sym.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP4, plot = F)
plot(Tol_Sym.lme_res_TP4)

##Model results
summary(Tol_Sym.lme_TP4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP4

REML criterion at convergence: -44.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5412 -0.7855 -0.1534  0.5518  2.2412 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.001295 0.03598 
 Residual             0.016289 0.12763 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.24607    0.04230  9.94436   5.818 0.000173 ***
OriginTransplant         0.04098    0.05210 42.00002   0.787 0.435975    
SiteSS                  -0.03092    0.05210 42.00002  -0.593 0.556032    
OriginTransplant:SiteSS -0.03256    0.07369 42.00002  -0.442 0.660816    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.616              
SiteSS      -0.616  0.500       
OrgnTrn:SSS  0.436 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP4)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.01 | [0.00, 1.00]
Site        |           0.04 | [0.00, 1.00]
Origin:Site |       4.63e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP4.res<-data.frame(summary(Tol_Sym.lme_TP4)$coefficients[-1,])
Tol_Sym_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP4.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP4)$Eta2)

Effect size of Origin for each Site

##KL
Tol_Sym.lme_TP4_KL<-lmer(log(Sym.prop+1)~Origin+(1|Genotype), data=Therm_TP4[which(Therm_TP4$Site=="KL"),])
summary(Tol_Sym.lme_TP4_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP4[which(Therm_TP4$Site == "KL"), ]

REML criterion at convergence: -26.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3745 -0.8456  0.0271  0.5432  2.0020 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002906 0.05391 
 Residual             0.012909 0.11362 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)       0.24607    0.04522  3.63655   5.442  0.00725 **
OriginTransplant  0.04098    0.04638 20.00000   0.884  0.38746   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.513

eta_squared(Tol_Sym.lme_TP4_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Sym.lme_TP4_SS<-lmer(log(Sym.prop+1)~Origin+(1|Genotype), data=Therm_TP4[which(Therm_TP4$Site=="SS"),])

boundary (singular) fit: see help('isSingular')

summary(Tol_Sym.lme_TP4_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP4[which(Therm_TP4$Site == "SS"), ]

REML criterion at convergence: -19.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.19030 -0.85642 -0.01412  0.53810  2.27217 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.00000  0.0000  
 Residual             0.01944  0.1394  
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.215142   0.040247 22.000000   5.346 2.29e-05 ***
OriginTransplant  0.008418   0.056918 22.000000   0.148    0.884    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Sym.lme_TP4_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |       9.93e-04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Sym_TP4.site.res<-data.frame(rbind(summary(Tol_Sym.lme_TP4_KL)$coefficients[-1,],
                                     summary(Tol_Sym.lme_TP4_SS)$coefficients[-1,]))
Tol_Sym_TP4.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Sym_TP4.site.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP4_KL)$Eta2, eta_squared(Tol_Sym.lme_TP4_SS)$Eta2)

##Combine results 
Tol_Sym_TP4.res<-rbind(Tol_Sym_TP4.res, Tol_Sym_TP4.site.res)
Tol_Sym_TP4.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP4.res))

##Summary statistics by Site and Origin
TP4_TolSym.sum<-summarySE(Therm_TP4, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP4_TolSym.plot<-ggplot(TP4_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1); TP4_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP4$Chl.prop)

shapiro.test(Therm_TP4$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP4$Chl.prop
W = 0.94204, p-value = 0.02131

#Not normal

hist(log(Therm_TP4$Chl.prop+1))

shapiro.test(log(Therm_TP4$Chl.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP4$Chl.prop + 1)
W = 0.9512, p-value = 0.04824

##Nearly normal

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP4<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP4)

boundary (singular) fit: see help('isSingular')

##Check residuals
Tol_Chl.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP4, plot = F)
plot(Tol_Chl.lme_res_TP4)

##Model results
summary(Tol_Chl.lme_TP4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP4

REML criterion at convergence: -139.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6397 -0.7647 -0.1492  0.5006  2.4251 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.000000 0.00000 
 Residual             0.001803 0.04246 
Number of obs: 47, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.09448    0.01226 43.00000   7.707 1.24e-09 ***
OriginTransplant        -0.01266    0.01734 43.00000  -0.730    0.469    
SiteSS                  -0.02659    0.01773 43.00000  -1.500    0.141    
OriginTransplant:SiteSS  0.02939    0.02479 43.00000   1.185    0.242    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.707              
SiteSS      -0.692  0.489       
OrgnTrn:SSS  0.494 -0.699 -0.715
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP4)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       6.26e-04 | [0.00, 1.00]
Site        |           0.02 | [0.00, 1.00]
Origin:Site |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP4.res<-data.frame(summary(Tol_Chl.lme_TP4)$coefficients[-1,])
Tol_Chl_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP4.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP4)$Eta2)

Effect size of Origin for each Site

##KL
Tol_Chl.lme_TP4_KL<-lmer(log(Chl.prop+1)~Origin+(1|Genotype), data=Therm_TP4[which(Therm_TP4$Site=="KL"),])

boundary (singular) fit: see help('isSingular')

summary(Tol_Chl.lme_TP4_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP4[which(Therm_TP4$Site == "KL"), ]

REML criterion at convergence: -76.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5133 -0.5558 -0.1258  0.5428  2.7208 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.000000 0.00000 
 Residual             0.001433 0.03785 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.09448    0.01093 22.00000   8.647 1.59e-08 ***
OriginTransplant -0.01266    0.01545 22.00000  -0.819    0.421    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP4_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_Chl.lme_TP4_SS<-lmer(log(Chl.prop+1)~Origin+(1|Genotype), data=Therm_TP4[which(Therm_TP4$Site=="SS"),])

boundary (singular) fit: see help('isSingular')

summary(Tol_Chl.lme_TP4_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin + (1 | Genotype)
   Data: Therm_TP4[which(Therm_TP4$Site == "SS"), ]

REML criterion at convergence: -64.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4873 -0.7211 -0.1933  0.5918  2.1830 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.000000 0.00000 
 Residual             0.002192 0.04681 
Number of obs: 23, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       0.06788    0.01412 21.00000   4.809 9.41e-05 ***
OriginTransplant  0.01673    0.01954 21.00000   0.856    0.402    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.722
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP4_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_Chl_TP4.site.res<-data.frame(rbind(summary(Tol_Chl.lme_TP4_KL)$coefficients[-1,],
                                     summary(Tol_Chl.lme_TP4_SS)$coefficients[-1,]))
Tol_Chl_TP4.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_Chl_TP4.site.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP4_KL)$Eta2, eta_squared(Tol_Chl.lme_TP4_SS)$Eta2)

##Combine results 
Tol_Chl_TP4.res<-rbind(Tol_Chl_TP4.res, Tol_Chl_TP4.site.res)
Tol_Chl_TP4.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP4.res))

##Summary statistics by Site and Origin
TP4_TolChl.sum<-summarySE(Therm_TP4, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP4_TolChl.plot<-ggplot(TP4_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 0.35); TP4_TolChl.plot

FvFm

#Check normality
hist(Therm_TP4$PAM.prop)

shapiro.test(Therm_TP4$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP4$PAM.prop
W = 0.86379, p-value = 5.159e-05

#Not normal

hist(log(Therm_TP4$PAM.prop+1))

shapiro.test(log(Therm_TP4$PAM.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP4$PAM.prop + 1)
W = 0.84113, p-value = 1.275e-05

##Still not normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP4<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP4)

##Check residuals
Tol_PAM.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP4, plot = F)
plot(Tol_PAM.lme_res_TP4)

##Model results
summary(Tol_PAM.lme_TP4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP4

REML criterion at convergence: -60.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1005 -0.4066 -0.0093  0.5802  2.4408 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.01070  0.1035  
 Residual             0.01028  0.1014  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)              0.74804    0.06651  2.73320  11.247  0.00228 **
OriginTransplant         0.03239    0.04138 42.00000   0.783  0.43818   
SiteSS                   0.02783    0.04138 42.00000   0.672  0.50502   
OriginTransplant:SiteSS -0.05803    0.05852 42.00000  -0.991  0.32713   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.311              
SiteSS      -0.311  0.500       
OrgnTrn:SSS  0.220 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP4)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       3.17e-04 | [0.00, 1.00]
Site        |       3.92e-05 | [0.00, 1.00]
Origin:Site |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP4.res<-data.frame(summary(Tol_PAM.lme_TP4)$coefficients[-1,])
Tol_PAM_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP4.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP4)$Eta2)

Effect size of Origin for each Site

##KL
Tol_PAM.lme_TP4_KL<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP4[which(Therm_TP4$Site=="KL"),])
summary(Tol_PAM.lme_TP4_KL)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP4[which(Therm_TP4$Site == "KL"), ]

REML criterion at convergence: -30.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.84954 -0.38862 -0.07926  0.54042  2.60201 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.011716 0.10824 
 Residual             0.009227 0.09606 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)       0.74804    0.06837  2.37263  10.941  0.00438 **
OriginTransplant  0.03239    0.03922 20.00001   0.826  0.41856   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.287

eta_squared(Tol_PAM.lme_TP4_KL)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##SS
Tol_PAM.lme_TP4_SS<-lmer(PAM.prop~Origin+(1|Genotype), data=Therm_TP4[which(Therm_TP4$Site=="SS"),])
summary(Tol_PAM.lme_TP4_SS)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin + (1 | Genotype)
   Data: Therm_TP4[which(Therm_TP4$Site == "SS"), ]

REML criterion at convergence: -25.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.08840 -0.38585  0.06028  0.59276  2.02689 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.008458 0.09197 
 Residual             0.012219 0.11054 
Number of obs: 24, groups:  Genotype, 3

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)       0.77587    0.06195  2.65242  12.524  0.00196 **
OriginTransplant -0.02563    0.04513 20.00000  -0.568  0.57634   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
OrgnTrnspln -0.364

eta_squared(Tol_PAM.lme_TP4_SS)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Origin    |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save results
Tol_PAM_TP4.site.res<-data.frame(rbind(summary(Tol_PAM.lme_TP4_KL)$coefficients[-1,],
                                     summary(Tol_PAM.lme_TP4_SS)$coefficients[-1,]))
Tol_PAM_TP4.site.res$Predictor<-c("KL Origin", "SS Origin")
Tol_PAM_TP4.site.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP4_KL)$Eta2, eta_squared(Tol_PAM.lme_TP4_SS)$Eta2)

##Combine results 
Tol_PAM_TP4.res<-rbind(Tol_PAM_TP4.res, Tol_PAM_TP4.site.res)
Tol_PAM_TP4.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP4.res))

##Summary statistics by Site and Origin
TP4_TolPAM.sum<-summarySE(Therm_TP4, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP4_TolPAM.plot<-ggplot(TP4_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.6, 1); TP4_TolPAM.plot

Save Results

##Combine Results
Tol_T4.res<-rbind(Tol_PAM_TP4.res, Tol_Sym_TP4.res, Tol_Chl_TP4.res)

##Add Timepoint
Tol_T4.res$TimeP<-rep("TP4", nrow(Tol_T4.res))

Effect Size over Time

##Dataframe of effect size results
Tol.ES<-rbind(Tol_T1.res, Tol_T2.res, Tol_T2.res, Tol_T4.res)
Tol.ES<-Tol.ES %>% dplyr::rename(p = Pr...t..)
Tol.ES<-Tol.ES[,c("TimeP", "Predictor", "Response", "EtaSq", "p")]

##Site specific results
Tol.ES<-Tol.ES[which(Tol.ES$Predictor=="KL Origin" | Tol.ES$Predictor=="SS Origin"),]
Tol.ES<-Tol.ES %>% separate_wider_delim(cols="Predictor", delim=" ", names=c("Site", "Predictor"), cols_remove = TRUE)

##Add Metric names
Tol.ES$Metric<-ifelse(Tol.ES$Response== "FvFm Retention", "FvFm", ifelse(
  Tol.ES$Response== "Symbiont Retention", "Symbionts", ifelse(
    Tol.ES$Response=="Chlorophyll Retention", "Chlorophyll", NA)))

##Add Significance levels
Tol.ES$Sig<-ifelse(Tol.ES$p<0.001, "***", ifelse(Tol.ES$p<0.01, "**", 
            ifelse(Tol.ES$p<0.05, "*", ifelse(Tol.ES$p<0.1, "-", NA))))

Symbionts

Tol_Sym.ES.plot<-ggplot(Tol.ES[which(Tol.ES$Metric=="Symbionts"),], aes(x=TimeP, y=EtaSq, fill=Site))+
  geom_bar(stat="identity", position=position_dodge())+
    scale_fill_manual(values=Site.colors.o)+
theme_classic()+
  ggtitle("Symbiont Retention")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
   labs(x="Time Point", y=expression(paste("Effect Size (p", eta^2, ")")), colour="Site")+
  ylim(0, 0.4)+
  geom_text(aes(label=Sig), vjust=-0.02, color="black", position=position_dodge(0.9), size=levels.sz, fontface="bold"); Tol_Sym.ES.plot

Chlorophyll

Tol_Chl.ES.plot<-ggplot(Tol.ES[which(Tol.ES$Metric=="Chlorophyll"),], aes(x=TimeP, y=EtaSq, fill=Site))+
  geom_bar(stat="identity", position=position_dodge())+
    scale_fill_manual(values=Site.colors.o)+
theme_classic()+
  ggtitle("Chlorophyll Retention")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
   labs(x="Time Point", y=expression(paste("Effect Size (p", eta^2, ")")), colour="Site")+
  ylim(0, 0.4)+
  geom_text(aes(label=Sig), vjust=-0.02, color="black", position=position_dodge(0.9), size=levels.sz, fontface="bold"); Tol_Chl.ES.plot

FvFm

Tol_PAM.ES.plot<-ggplot(Tol.ES[which(Tol.ES$Metric=="FvFm"),], aes(x=TimeP, y=EtaSq, fill=Site))+
  geom_bar(stat="identity", position=position_dodge())+
    scale_fill_manual(values=Site.colors.o)+
theme_classic()+
  ggtitle("Photochemical Efficiency Retention")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
   labs(x="Time Point", y=expression(paste("Effect Size (p", eta^2, ")")), colour="Site")+
  ylim(0, 0.4)+
  geom_text(aes(label=Sig), vjust=-0.02, color="black", position=position_dodge(0.9), size=levels.sz, fontface="bold"); Tol_PAM.ES.plot

Trade Offs

Initial T1 and T2

Average Performance

Calculating the average of Growth (TP1 to TP2) and Thermal Tolerance (Retention of Symbionts, Chlorophyll, and Photochemical Efficiency at TP1) for each Set (Site, Genotype, Origin).

Growth

##Average
names(Growth_T1.2)

 [1] "Site"       "ID"         "Orig"       "Origin"     "Genotype"   "TL.TP1_cm"  "TL.TP2_cm" 
 [8] "TL.TP3_cm"  "TL.TP4_cm"  "Set"        "Site.Orig"  "TP1"        "TP2"        "TP1.2_days"
[15] "Ext_cm"     "TLE_cm.day"

T1.2_Growth_Set<-aggregate(Growth_T1.2$TLE_cm.day, list(Growth_T1.2$Site, Growth_T1.2$Genotype, Growth_T1.2$Orig), mean)
names(T1.2_Growth_Set)<-c("Site", "Genotype", "Orig", "TLE_cm.day")

##Add Standard Error
T1.2_Growth_Set$TLE_se<-aggregate(Growth_T1.2$TLE_cm.day, list(Growth_T1.2$Site, Growth_T1.2$Genotype, Growth_T1.2$Orig), std.error)[,4]

Tolerance

##Average
names(Therm_TP1)

 [1] "TimeP"         "Site"          "Genotype"      "Origin"        "ID"            "RandN"        
 [7] "Treat"         "Treatment"     "Orig"          "Set"           "Site.Orig"     "SA_cm2"       
[13] "Chl_ug.cm2"    "Sym10.6_cm2"   "Fv_Fm"         "Chl_ug.cm2_C"  "Sym10.6_cm2_C" "Fv_Fm_C"      
[19] "Chl.prop"      "Sym.prop"      "PAM.prop"     

T1_Therm_Set<-aggregate(Therm_TP1[,c(19:21)], list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), mean)
names(T1_Therm_Set)<-c("Site", "Genotype", "Orig", "Chl.prop", "Sym.prop", "PAM.prop")

##Add Standard Error
T1_Therm_Set$Chl_se<-aggregate(Therm_TP1$Chl.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]

T1_Therm_Set$Sym_se<-aggregate(Therm_TP1$Sym.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]

T1_Therm_Set$PAM_se<-aggregate(Therm_TP1$PAM.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]

Merge Growth and Tolerance

T1.2_Perf_Set<-merge(T1.2_Growth_Set, T1_Therm_Set, all=TRUE)
T1.2_Perf_Set$Set<-paste(T1.2_Perf_Set$Site, T1.2_Perf_Set$Genotype, T1.2_Perf_Set$Orig, sep=".")

##Set factor variables
T1.2_Perf_Set$Site<-factor(T1.2_Perf_Set$Site, levels=c("KL", "SS"))
T1.2_Perf_Set$Genotype<-factor(T1.2_Perf_Set$Genotype, levels=c("AC8", "AC10", "AC12"))
T1.2_Perf_Set$Orig<-factor(T1.2_Perf_Set$Orig, levels=c("N", "T"))

##Add a Sample Set Variable
T1.2_Perf_Set$Set<-paste(T1.2_Perf_Set$Site, T1.2_Perf_Set$Genotype, T1.2_Perf_Set$Orig, sep=".")
T1.2_Perf_Set$Set<-factor(T1.2_Perf_Set$Set)

##Add Site.Orig variable
T1.2_Perf_Set$Site.Orig<-paste(T1.2_Perf_Set$Site, T1.2_Perf_Set$Orig, sep=".")
T1.2_Perf_Set$Site.Orig<-factor(T1.2_Perf_Set$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))

Correlation between Growth and Tolerance

Symbiont Retention

cor.test(T1.2_Perf_Set$TLE_cm.day, T1.2_Perf_Set$Sym.prop, method="spearman")

    Spearman's rank correlation rho

data:  T1.2_Perf_Set$TLE_cm.day and T1.2_Perf_Set$Sym.prop
S = 462, p-value = 0.03733
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.6153846 

T1.2_TLE_Sym.plot<-ggplot(T1.2_Perf_Set, aes(x=TLE_cm.day, y=Sym.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Symbiont Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = -.15, label=expression(bolditalic(paste(r[S], " = -0.615, p = 0.037"))), size=sig.sz, hjust = 0); T1.2_TLE_Sym.plot

Chlorophyll Retention

cor.test(T1.2_Perf_Set$TLE_cm.day, T1.2_Perf_Set$Chl.prop, method="spearman")

    Spearman's rank correlation rho

data:  T1.2_Perf_Set$TLE_cm.day and T1.2_Perf_Set$Chl.prop
S = 454, p-value = 0.04884
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.5874126 

T1.2_TLE_Chl.plot<-ggplot(T1.2_Perf_Set, aes(x=TLE_cm.day, y=Chl.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Chlorophyll Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = -.05, label=expression(bolditalic(paste(r[S], " = -0.587, p = 0.049"))), size=sig.sz, hjust = 0); T1.2_TLE_Chl.plot

FvFm Retention

cor.test(T1.2_Perf_Set$TLE_cm.day, T1.2_Perf_Set$PAM.prop, method="spearman")

    Spearman's rank correlation rho

data:  T1.2_Perf_Set$TLE_cm.day and T1.2_Perf_Set$PAM.prop
S = 480, p-value = 0.01883
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.6783217 

T1.2_TLE_PAM.plot<-ggplot(T1.2_Perf_Set, aes(x=TLE_cm.day, y=PAM.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Photo. Effic. Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.7, label=expression(bolditalic(paste(r[S], " = -0.678, p = 0.019"))), size=sig.sz, hjust = 0); T1.2_TLE_PAM.plot

Final T3 and T4

Average Performance

Calculating the average of Growth (TP3 to TP4) and Thermal Tolerance (Retention of Symbionts, Chlorophyll, and Photochemical Efficiency at TP3) for each Set (Site, Genotype, Origin).

Growth

##Average
names(Growth_T3.4)

 [1] "Site"       "ID"         "Orig"       "Origin"     "Genotype"   "TL.TP1_cm"  "TL.TP2_cm" 
 [8] "TL.TP3_cm"  "TL.TP4_cm"  "Set"        "Site.Orig"  "TP3"        "TP4"        "TP3.4_days"
[15] "Ext_cm"     "TLE_cm.day"

T3.4_Growth_Set<-aggregate(Growth_T3.4$TLE_cm.day, list(Growth_T3.4$Site, Growth_T3.4$Genotype, Growth_T3.4$Orig), mean)
names(T3.4_Growth_Set)<-c("Site", "Genotype", "Orig", "TLE_cm.day")

##Add Standard Error
T3.4_Growth_Set$TLE_se<-aggregate(Growth_T3.4$TLE_cm.day, list(Growth_T3.4$Site, Growth_T3.4$Genotype, Growth_T3.4$Orig), std.error)[,4]

Tolerance

##Average
names(Therm_TP4)

 [1] "TimeP"         "Site"          "Genotype"      "Origin"        "ID"            "RandN"        
 [7] "Treat"         "Treatment"     "Orig"          "Set"           "Site.Orig"     "SA_cm2"       
[13] "Chl_ug.cm2"    "Sym10.6_cm2"   "Fv_Fm"         "Chl_ug.cm2_C"  "Sym10.6_cm2_C" "Fv_Fm_C"      
[19] "Chl.prop"      "Sym.prop"      "PAM.prop"     

T4_Therm_Set<-aggregate(Therm_TP4[,c(19:21)], list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), mean)
names(T4_Therm_Set)<-c("Site", "Genotype", "Orig", "Chl.prop", "Sym.prop", "PAM.prop")

##Add Standard Error
T4_Therm_Set$Chl_se<-aggregate(Therm_TP4$Chl.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]

T4_Therm_Set$Sym_se<-aggregate(Therm_TP4$Sym.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]

T4_Therm_Set$PAM_se<-aggregate(Therm_TP4$PAM.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]

Merge Growth and Tolerance

T3.4_Perf_Set<-merge(T3.4_Growth_Set, T4_Therm_Set, all=TRUE)
T3.4_Perf_Set$Set<-paste(T3.4_Perf_Set$Site, T3.4_Perf_Set$Genotype, T3.4_Perf_Set$Orig, sep=".")

##Set factor variables
T3.4_Perf_Set$Site<-factor(T3.4_Perf_Set$Site, levels=c("KL", "SS"))
T3.4_Perf_Set$Genotype<-factor(T3.4_Perf_Set$Genotype, levels=c("AC8", "AC10", "AC12"))
T3.4_Perf_Set$Orig<-factor(T3.4_Perf_Set$Orig, levels=c("N", "T"))

##Add a Sample Set Variable
T3.4_Perf_Set$Set<-paste(T3.4_Perf_Set$Site, T3.4_Perf_Set$Genotype, T3.4_Perf_Set$Orig, sep=".")
T3.4_Perf_Set$Set<-factor(T3.4_Perf_Set$Set)

##Add Site.Orig variable
T3.4_Perf_Set$Site.Orig<-paste(T3.4_Perf_Set$Site, T3.4_Perf_Set$Orig, sep=".")
T3.4_Perf_Set$Site.Orig<-factor(T3.4_Perf_Set$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))

Correlation between Growth and Tolerance

Symbiont Retention

cor.test(T3.4_Perf_Set$TLE_cm.day, T3.4_Perf_Set$Sym.prop, method="spearman")

    Spearman's rank correlation rho

data:  T3.4_Perf_Set$TLE_cm.day and T3.4_Perf_Set$Sym.prop
S = 444, p-value = 0.06663
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.5524476 

T3.4_TLE_Sym.plot<-ggplot(T3.4_Perf_Set, aes(x=TLE_cm.day, y=Sym.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Symbiont Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.1, label=expression(italic(paste(r[S], " = -0.552, p = 0.067"))), size=sig.sz, hjust = 0); T3.4_TLE_Sym.plot

Chlorophyll Retention

cor.test(T3.4_Perf_Set$TLE_cm.day, T3.4_Perf_Set$Chl.prop, method="spearman")

    Spearman's rank correlation rho

data:  T3.4_Perf_Set$TLE_cm.day and T3.4_Perf_Set$Chl.prop
S = 148, p-value = 0.327
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.3272727 

T3.4_TLE_Chl.plot<-ggplot(T3.4_Perf_Set, aes(x=TLE_cm.day, y=Chl.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Chlorophyll Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.065, label=expression(italic(paste(r[S], " = 0.327, p = 0.327"))), size=sig.sz, hjust = 0); T3.4_TLE_Chl.plot

FvFm Retention

cor.test(T3.4_Perf_Set$TLE_cm.day, T3.4_Perf_Set$PAM.prop, method="spearman")

    Spearman's rank correlation rho

data:  T3.4_Perf_Set$TLE_cm.day and T3.4_Perf_Set$PAM.prop
S = 500, p-value = 0.007353
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.7482517 

T3.4_TLE_PAM.plot<-ggplot(T3.4_Perf_Set, aes(x=TLE_cm.day, y=PAM.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Photo. Effic. Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.5, label=expression(bolditalic(paste(r[S], " = -0.748, p = 0.007"))), size=sig.sz, hjust = 0); T3.4_TLE_PAM.plot

Figures

Figure 4 Growth

Adjust Figures for Panel

##TP 1 to TP 2
TP1.2_Growth.plot<-TP1.2_Growth.plot+ ggtitle("TP1 to TP2")+
  theme(plot.title = element_text(colour="black", 
                                  size=panel.lab.sz, hjust=0.5), 
        legend.position=c(0.25, 0.85), legend.direction="horizontal")+
  labs(colour=NULL)+
  guides(color=guide_legend(nrow=2, byrow=FALSE)) 


##TP 3 to TP 4
TP3.4_Growth.plot<-TP3.4_Growth.plot+ ggtitle("TP3 to TP4")+
  theme(plot.title = element_text(colour="black", 
                                  size=panel.lab.sz, hjust=0.5), 
        legend.position="none")

Figure 4

##Create Panel
Growth_fig<-plot_grid(TP1.2_Growth.plot, TP3.4_Growth.plot, 
                    rel_widths=1, rel_heights = 1,
                    nrow=1, ncol=2, byrow=T, labels = c("A", "B"), 
                    label_size=panel.lab.sz)

##Save Figure
ggsave(filename="Figures/03_Performance/Fig4_Growth.png", plot=Growth_fig, dpi=300, width=10, height=4.5, units="in")

Figure 5 Thermal Tolerance

Adjust Figures for Panel

##FvFm Retention
TP1_TolPAM.plot<-TP1_TolPAM.plot + ggtitle("TP1")+ 
  labs(x="", y="Photo. Effic. Retention", colour=NULL)+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position=c(0.5, 0.15), legend.direction="horizontal")+
  guides(color=guide_legend(nrow=2, byrow=FALSE))

TP2_TolPAM.plot<-TP2_TolPAM.plot+ ggtitle("TP2")+ labs(x="", y="") +
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position="none")

TP3_TolPAM.plot<-TP3_TolPAM.plot+ ggtitle("TP3")+ labs(x="", y="") +
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position="none")

TP4_TolPAM.plot<-TP4_TolPAM.plot+ ggtitle("TP4")+ labs(x="", y="") +
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position="none")

Tol_PAM.ES.plot<-Tol_PAM.ES.plot + ggtitle("Photo. Effic.") + labs(x="") +
  theme(legend.position=c(.8, .8))


##Symbiont Retention
TP1_TolSym.plot<- TP1_TolSym.plot + labs(x="") + theme(legend.position="none")

TP2_TolSym.plot<- TP2_TolSym.plot + labs(x="", y= "") + theme(legend.position="none")

TP3_TolSym.plot<- TP3_TolSym.plot + labs(x="", y= "") + theme(legend.position="none")

TP4_TolSym.plot<- TP4_TolSym.plot + labs(x="", y= "") + theme(legend.position="none")

Tol_Sym.ES.plot<-Tol_Sym.ES.plot +  ggtitle("Symbionts") + labs(x="") + theme(legend.position="none")


##Chlorophyll Retention
TP1_TolChl.plot<-TP1_TolChl.plot + theme(legend.position="none")

TP2_TolChl.plot<-TP2_TolChl.plot + labs(y="") + theme(legend.position="none")

TP3_TolChl.plot<-TP3_TolChl.plot + labs(y="") + theme(legend.position="none")

TP4_TolChl.plot<-TP4_TolChl.plot + labs(y="") + theme(legend.position="none")

Tol_Chl.ES.plot<-Tol_Chl.ES.plot + ggtitle("Chlorophyll")+ theme(legend.position="none")

Figure 5

##Create Panel
Tolerance_fig<-plot_grid(TP1_TolPAM.plot, TP2_TolPAM.plot, 
                      TP3_TolPAM.plot, TP4_TolPAM.plot, Tol_PAM.ES.plot, 
                      TP1_TolSym.plot, TP2_TolSym.plot, 
                      TP3_TolSym.plot, TP4_TolSym.plot, Tol_Sym.ES.plot, 
                      TP1_TolChl.plot, TP2_TolChl.plot, 
                      TP3_TolChl.plot, TP4_TolChl.plot, Tol_Chl.ES.plot, 
                    rel_widths=c(0.8, 0.8, 0.8, 0.8, 1), rel_heights = 1,
                    nrow=3, ncol=5, byrow=T, labels = NULL)

Warning: Removed 5 rows containing missing values (`geom_text()`).Warning: Removed 7 rows containing missing values (`geom_text()`).Warning: Removed 8 rows containing missing values (`geom_text()`).

##Save Figure
ggsave(filename="Figures/03_Performance/Fig5_Tolerance.png", plot=Tolerance_fig, dpi=300, width=16, height=12, units="in")

Figure 6 Growth vs Tolerance

Adjust Figures for Panel

##Timepoint 1 to 2
T1.2_TLE_PAM.plot<-T1.2_TLE_PAM.plot + ggtitle("Photo. Effic.")+ labs(x="", y = "Tolerance")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), legend.position="none")

T1.2_TLE_Sym.plot<- T1.2_TLE_Sym.plot + ggtitle("Symbionts")+ labs(x="", y = "")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), legend.position="none") + ylim(-.15, .8)

T1.2_TLE_Chl.plot<- T1.2_TLE_Chl.plot + ggtitle("Chlorophyll")+ labs(x="", y = "")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), legend.position="right")

##Timepoint 3 to 4
T3.4_TLE_PAM.plot<-T3.4_TLE_PAM.plot + labs(x="Growth", y = "Tolerance") +
  theme(legend.position="none")

T3.4_TLE_Sym.plot<-T3.4_TLE_Sym.plot + labs(x="Growth", y = "") +
  theme(legend.position="none")

T3.4_TLE_Chl.plot<-T3.4_TLE_Chl.plot + labs(x="Growth", y = "") +
  theme(legend.position="right") 

Figure 6

##Create Panel
TradeOff_fig<-plot_grid(T1.2_TLE_PAM.plot, T1.2_TLE_Sym.plot, T1.2_TLE_Chl.plot,
                      T3.4_TLE_PAM.plot, T3.4_TLE_Sym.plot, T3.4_TLE_Chl.plot,
                    rel_widths=c(.7, .7, 1, .7, .7, 1), rel_heights = c(1, .9),
                    nrow=2, ncol=3, byrow=T, 
                    labels = c("A", "B", "C", "D", "E", "F"), 
                    label_size=panel.lab.sz)

`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).Warning: Removed 1 rows containing missing values (`geom_point()`).Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).Warning: Removed 1 rows containing missing values (`geom_point()`).Warning: is.na() applied to non-(list or vector) of type 'expression'

##Save Figure
ggsave(filename="Figures/03_Performance/Fig6_TradeOff.png", plot=TradeOff_fig, dpi=300, width=12, height=8, units="in")

Figure S4 Initial Heat Assay

Adjust Figures for Panel

##FvFm
IN_PAM.plot<-IN_PAM.plot + labs(x="", y="Fv/Fm") + theme(legend.position="none") + ylim(0.5, 0.655)

IN_TolPAM.plot<-IN_TolPAM.plot + labs(x="", y="Retention") + 
  theme(legend.position="none") + ylim(0.75, 1)

##Symbionts
IN_Sym.plot<-IN_Sym.plot + labs(x="", y="Symbionts") + theme(legend.position="none") + ylim(0,1)

IN_TolSym.plot<-IN_TolSym.plot + labs(x="", y="Retention") + 
  theme(legend.position="none") + ylim(.25, .75)

##Chlorophyll
IN_Chl.plot<-IN_Chl.plot + labs(x="Site and Treatment", y="Chlorophyll") + 
  theme(legend.position="none") + ylim(0,1.25)

IN_TolChl.plot<-IN_TolChl.plot + labs(x="Site", y="Retention") + 
  theme(legend.position="none") + ylim(0, 0.25)

Figure S4

##Create Panel
IN_HeatAssay_fig<-plot_grid(IN_PAM.plot, IN_TolPAM.plot, 
                      IN_Sym.plot, IN_TolSym.plot,
                      IN_Chl.plot, IN_TolChl.plot,
                    rel_widths=c(1, 0.75), 
                    rel_heights = c(1, 1, 1),
                    nrow=3, ncol=2, byrow=T, 
                    labels = c("A", "B", "C", "D", "E", "F"), 
                    label_size=panel.lab.sz)

##Save Figure
ggsave(filename="Figures/03_Performance/FigS4_InitialHeatAssay.png", plot=IN_HeatAssay_fig, dpi=300, width=8, height=11, units="in")

Tables

Table S5 Performance LM

##Combine Results Tables
TableS5A_Growth<-data.frame(rbind(TLE_T1.2.res, TLE_T3.4.res))
TableS5B_Tolerance<-data.frame(rbind(Tol_T1.res, Tol_T2.res, Tol_T2.res, Tol_T4.res))
names(TableS5B_Tolerance)[9]<-"Timepoint"

TableS5_Perf.LM<-rbind(TableS5A_Growth, TableS5B_Tolerance)
  
##Organize
names(TableS5_Perf.LM)

[1] "Estimate"   "Std..Error" "df"         "t.value"    "Pr...t.."   "Predictor" 
[7] "EtaSq"      "Response"   "Timepoint" 

TableS5_Perf.LM<-TableS5_Perf.LM %>% dplyr::rename(SE = Std..Error, DF = df, t = t.value, p = Pr...t.., EffectSize = EtaSq)
TableS5_Perf.LM<-TableS5_Perf.LM[,c("Timepoint", "Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EffectSize")]

#Round to 3 digits
TableS5_Perf.LM$Estimate<-round(TableS5_Perf.LM$Estimate, 3)
TableS5_Perf.LM$SE<-round(TableS5_Perf.LM$SE, 3)
TableS5_Perf.LM$t<-round(TableS5_Perf.LM$t, 3)
TableS5_Perf.LM$p<-round(TableS5_Perf.LM$p, 3)
TableS5_Perf.LM$EffectSize<-round(TableS5_Perf.LM$EffectSize, 3)

#Integer
TableS5_Perf.LM$DF<-round(TableS5_Perf.LM$DF, 0)

##Write Out Table
write.csv(TableS5_Perf.LM, "Tables/TableS5_Performance_LM_Results.csv", row.names=FALSE)

Table S6 Initial Heat Assay LM

##Combine Results Tables
TableS6_IN_Tolerance<-data.frame(rbind(PAM_IN.res, Sym_IN.res, Chl_IN.res))

##Organize
names(TableS6_IN_Tolerance)

[1] "estimate"  "SE"        "df"        "t.ratio"   "p.value"   "Predictor"
[7] "Response"  "EtaSq"    

TableS6_IN_Tolerance<-TableS6_IN_Tolerance %>% dplyr::rename(Estimate = estimate, DF = df, t = t.ratio, p = p.value, EffectSize = EtaSq)
TableS6_IN_Tolerance<-TableS6_IN_Tolerance[,c("Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EffectSize")]

#Round to 3 digits
TableS6_IN_Tolerance$Estimate<-round(TableS6_IN_Tolerance$Estimate, 3)
TableS6_IN_Tolerance$SE<-round(TableS6_IN_Tolerance$SE, 3)
TableS6_IN_Tolerance$t<-round(TableS6_IN_Tolerance$t, 3)
TableS6_IN_Tolerance$p<-round(TableS6_IN_Tolerance$p, 3)
TableS6_IN_Tolerance$EffectSize<-round(TableS6_IN_Tolerance$EffectSize, 3)

#Integer
TableS6_IN_Tolerance$DF<-round(TableS6_IN_Tolerance$DF, 0)

##Write Out Table
write.csv(TableS6_IN_Tolerance, "Tables/TableS6_Initial_Tolerance_LM_Results.csv", row.names=FALSE)